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Abstract 


Parameter sensitivity is defined as the estimation of changes in the modeling 
functions and design variables due to small changes in the fixed parameters of the 
formulation. There are currently several methods for estimating parameter sensitivities 
which either require difficult to obtain second order information, or do not return reliable 
estimates for the derivatives. Additionally, all the methods assume that the set of active 
constraints does not change in a neighborhood of the estimation point. If the active set 
does in fact change, than any extrapolations based on these derivatives may be in error. It 
is the objective of this work to investigate new methods for estimating parameter 
sensitivities that are more efficient than current methods for estimating sensitivities when 
the active set changes. 

The new method proposed for estimating sensitivity derivatives is based on the 
recursive quadratic programming (RQP) method and in conjunction a differencing formula 
to produce estimates of the sensitivities. This method is compared to existing methods and 
is shown to be very competitive in terms of the number of function evaluations required. 

In terms of accuracy, the method is shown to be equivalent to a modified version of the 
Kuhn-Tucker method, where the Hessian of the Lagrangian is estimated using the BFS 
method employed by the RQP algorithm. Initial testing on a test set with known 
sensitivities demonstrates that the method can accurately calculate the parameter sensitivity. 

To handle changes in the active set, a deflection algorithm is proposed for those 
cases where the new actives set of constraints remains linearly independent. For those 
cases where dependencies occur, a directional derivative is proposed. A few simple 
examples are included for the algorithm, but extensive testing has not yet been performed. 
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1 . 


Introduction 


Estimation of the sensitivity of problem functions with respect to problem variables 
forms the basis for many of our modem day algorithms for engineering optimization. The 
most common application of problem sensitivities has been in the calculation of objective 
function and constraint partial derivatives for determining search directions and optimality 
conditions. A second form of sensitivity analysis, parameter sensitivity, has also become 
an important topic in recent years with the advent of renewed research in the optimization of 
large engineering systems by means of decomposition methods. By parameter sensitivity, 
we refer to the estimation of changes in the modeling functions and current design variables 
due to small changes in the fixed parameters of the formulation. Methods for calculating 
these derivatives have been proposed and have been used as the basis of a method for 
multi-level decomposition of large engineering problems [Sobieski, 1982]. Two 
drawbacks to estimating parameter sensitivities by current methods have been: (1) the need 
for second order information about the Lagrangian at the cuiTent point, and (2) the 
estimates assume no change in the active set of constraints. The objectives of this work 
were to investigate solutions to these two problems. 


1.1. STANDARD NOTATION 

To provide a framework about which we can discuss the various ways sensitivity 
analysis can be performed, the following standard form of the nonlinear programming 
problem, which explicitly represents the problem parameters, is presented. 


Minimize: f(x,P) 

Objective function 

(1.1) 

Subject to: hi(x,P) = 0 

Equality constraints 1=1,L 

(1.2) 

gj(x,P) > 0 

Inequality constraints j = 1,J 

(1.3) 

x min < x < x max 

Variable bounds 

(1.4) 

x = (xi,X2,...,x n ) 

Design variables 

(1.5) 

P = (Pl>P2v.,Pk) 

Problem parameters 

(1-6) 


In the above formulation, we assume that the problem functions f, g, and h can be 
either linear or nonlinear functions of the design variables. We also assume that the 
problem parameters P, are held fixed during the course of the optimization. Any candidate 
solution point, x* , must satisfy the following first order Kuhn-Tucker conditions: 


V x L(x,v,u) = 0 


(1.7) 

hi(x) = 0 

1= 1,L 

(1.8) 

gj(x) > 0 

j = U 

(1.9) 


( 1 . 10 ) 

( 1 . 11 ) 

( 1 . 12 ) 

At some point, usually the optimal point, we are interested in understanding the 
effect that changes in P will have on our proposed solution x*. Therefore we seek the 
sensitivities, df/dP, dx/dP, and d(h,g)/3P 1 . In this report, we will propose a new 
algorithm based on the Recursive Quadratic Programming (RQP) method for estimating 
these parameter sensitivities. The following sections provide a description of this algorithm 
and how it relates to current methods, a discussion of the implementation issues, and some 
initial testing on a test set of known characteristics. In addition, section 6 proposes some 
solutions for estimating sensitivities in those cases where the active set of the constraints 
changes when the parameter is changed. 

1 The notation (h,g) refers to the set of constraints active at the current point 


ujgj(x) = 0 j = U 
uj > 0 j = l,J 

where the Lagrangian L, is given by: 

L(x,v,u) = f(x) + Zvi hi(x) - Xu; gj(x) 
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2 . Background 


The standard problem of parameter sensitivity analysis is to indicate how the 
objective function, constraints, and optimum design variables will change when problem 
parameters or design variables are changed from their current values. Parameter Sensitivity 
analysis is usually performed at a candidate optimum point where we might be interested in 
studying how the optimal design might be effected by changes in specifications, variability 
due to manufacturing, or operational noises. In this chapter we present a historical 
overview of the significant developments in sensitivity analysis and provide a review and 
assessment of current parameter sensitivity methods. The final section of the chapter 
reviews work done in estimating parameter sensitivities for those cases where the active 
constraint set changes. 

2. 1 . REVIEW OF PARAMETER SENSITIVITY METHODS 


The roots of sensitivity analysis can be traced to Lagrange (1881) when he 
suggested solving equality constrained extrema problems by finding the solution x*,and 
v*, for the equations 

V x L(x,v) = 0 (2.1) 

h(x) = 0 (2.2) 

where 

L(x,v) = f(x) + Zvjh](x) (2.3) 

where the vi are undetermined multipliers or Lagrange multipliers. The paper did not 
provide the conditions for when solutions of equation (2. 1-2.3), were actual solutions of 
the extrema problems or how to interpret the Lagrange multipliers. 

Samuelson (1947) gave several interpretations of Lagrange multipliers in an 
economic setting. He developed approaches based on using Lagrange multipliers to solve 
different economic models and was the first to clearly identify Lagrange multipliers as 
shadow prices in an economic context. Kuhn and Tucker (1951) presented conditions for 
relative extrema which use the Lagrange multipliers to establish optimality (ref. eq. 1.7 - 
1.12). Since 1951 several constraint qualifications and extensions to these conditions have 
been proposed and are described in Bazaraa and Shetty (1979). 

Dantzig (1963) brought forth the idea of "Post Optimality Analysis" for linear 

programs. Dantzig described post optimality analysis as the calculation of the sensitivity of 

the optimum with respect to changes in the problem parameters. Sensitivity analysis has 

been widely used in linear programming, a good survey of its use is provided by Gal 
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(1984). 


Fiacco et al. (1968,1974,1976,1983) has also done extensive research in the area of 
sensitivity analysis. His book "Introduction to Sensitivity and Stability Analysis" (1983) 
covers the significant developments in the field of sensitivity analysis prior to 1982. He 
has published many articles on sensitivity analysis, and has probably been the most active 
researcher of sensitivity analysis for nonlinear programming problems. 

In the following subsections, we will discuss past work related to the determination 
of sensitivity information for nonlinear programming problems. The methods we will 
discuss range from the most simplistic approach of reoptimization to more elaborate 
approaches based on the Kuhn-Tucker conditions or advanced optimization methods. 

2.1.1. Brute Force Methods 

The simplest, and probably most used method, for parameter sensitivity analysis is 
to re-optimize the problem for the new values of the problem parameters and plot the 
trends. We will refer to this as the Brute Force method. The Brute Force method is 
probably the most accurate of the methods available (for large variations in Ap, but can 
experience round off and truncation errors when used to approximate derivatives) but it can 
be computationally expensive even for small problems. Examples of its use in the literature 
are given in Arbuckle and Sliwa (1984) and Robertson and Gabriele (1987). 

Armacost and Fiacco (1974) and McKeown (1980 b) describe a direct approach to 
calculating parameter sensitivities based on the central difference approximation given 
below 


df* f(x*,p + Ap) - f(x*,p - Ap) 

dP = 2Ap (2 ‘ 4) 

ax* x*(p + A p) - x*(p - Ap) 5) 

2Ap 

This method requires the problem to be reoptimized (to a high degree of accuracy) for two 
different values of the parameter. McKeown states that this method should not be used as a 
primary method for the calculation of sensitivities because it is computationally expensive. 

2.1.2. Kuhn-Tucker Methods 

To avoid the computational expense of reoptimization, several researchers have 
developed sensitivity methods based on the Kuhn-Tucker conditions (1.7) - (1.12). Two 
types of algorithms have resulted, those that differentiate the Kuhn-Tucker conditions with 
respect to p, and those that differentiate the optimality conditions for penalty functions. 
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In the former category, a set of Kuhn-Tucker sensitivity equations have been 
derived independently by several authors (Armacost and Fiacco 1974, Sobieski et. al. 
1981, McKeown 1980 b) and result in the following linear system of equations. 


V* L -V x (h,g) 

dx 


rdv xL i 

X v 

dpi 


dpi 

_ V x (h,g)T 0 _ 

d(v,u) 


d(h,g) 

- dpi - 


- dpi _ 


This linear system can be solved for the sensitivity of the design variables with 
respect to a problem parameter dx/dpj, and the sensitivity of the Lagrange multipliers with 
respect to pi, d(v,u)/dpj. These can then be used to determine the sensitivity of the 
objective function with respect to pi by the following 


df__ df af T ax 
dpi 5pi Sx" <3pT 


(2.7) 


For any change in the parameter Api, the new optimum value of the objective 
function or design variables can be estimated from the linear extrapolations 


fnew - f( x *old) + Apj 

* * . dx 

x new - x old + Api ^7 


(2.8) 

(2.9) 


These equations are bounded by the assumption that the active set remains the 
same. An estimate of when the active set will change can be made by examining the 
Lagrange multipliers of the active inequality constraints and linear approximations of the 
inactive constraints. An inequality constraint should leave the active set when its Lagrange 
multiplier goes to zero. The corresponding value of Api where this occurs is predicted by 
using the linear prediction 

ui 

Apj = — J — j e active set of constraints (2.10) 

v a PL 

A new inequality constraint will enter the active set when its value goes to zero. A linear 
prediction for when this happens is given by 


Api = 



£i 

dgj T dx^j 
dx ^PiJ 


4 active set of constraints 


( 2 . 11 ) 
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We can predict the change in active set to occur at the smallest value of Apj obtained from 
applying equations 2.10 and 2.1 1 to all constraints. 

Fiacco (1974,1980,1983) has developed first and second order extrapolation 
techniques to predict the new value of the optimum when parameters are perturbed. 
Armacost and Fiacco have developed a second order extrapolation for the objective function 
value for the special case where the problem parameters are confined to being the right hand 
side values of the constraints. This provides second order response information for the 
objective function using the Lagrange multipliers and the parti als with respect to P of the 
Lagrange multipliers. 

Sobieski, et. al. (1981) observed that a more accurate estimate of f n ew given in 
(2.8) can be obtained if the value of x new given in (2.9) is used to calculate the value of the 
objective function at a perturbation Ap,. This will be a more accurate estimate for problems 
where the constraints are well behaved and not highly nonlinear, but the objective function 
is nonlinear. 

Barthelemy and Sobieski (1983) derived the following formula that can also be 
used to calculate the sensitivity of the objective function without the need to calculate 
9x*/3p, 

nine q 

+ < 2 - 12 > 
j=l 


The formula can be derived by assuming that objective function behaves like the 
Lagrangian in the region of the optimum. This formula has also been derived by Fiacco 
(1983) and McKeown (1980 b). 

Diewart (1984) has developed some new sensitivity theories for dealing with the 
addition of constraints at the solution of economic models before the solution of the 
sensitivity equations. This analysis is important because there may be short term 
restrictions on modifications that can be made to the system. The paper presents a 
recursive relationship that can be used to avoid refactorizing the sensitivity equations when 
a new constraint is added to the problem. The paper also presents equations that can be 
used to calculate a second order estimate of the location of the optimum, but this formula 
requires third order derivatives which are seldom available in engineering. 

2.1.3. Methods Based on the Extended Design Space 

Vanderplaats (1984 a, 1984 b) and Vanderplaats and Yoshida (1985, 1986) have 
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developed an approach for calculating the sensitivity based on the method of feasible 
directions. The sensitivities are estimated by extending the set of design variables to 
include the problem parameters for which a feasible direction is then determined. This 
method is known as the Extended Design Space (EDS) method. Of the methods discussed, 
it has the dual advantages of simplicity and efficiency. Vanderplaats (1984 a) reports that 
the EDS method can handle near active constraints, and is able to leave constraint 
linearizations. However, the method does suffer from a sensitivity to one of its algorithm 
parameters as reported in Vanderplaats and Cai (1987), and is unable to predict when 
constraints will leave the active set. The EDS method is also sensitive to the restriction of 
the move vector to be of length one. 

The EDS method can be used to assess the effect of perturbing several parameters at 
the same time. It is also able to solve for sensitivities of degenerate optimal points where 
either strict complementarity does not hold, or the constraint gradients of the active 
constraints are linearly dependent. The method seems to give good estimations for medium 
sized perturbations of the parameters, but for small perturbations the the Kuhn-Tucker 
method described above gives better results. Vanderplaats and Cai (1987) also report that 
there are some cases where the EDS algorithm can produce incorrect values of the 
sensitivity derivatives. 

Vanderplaats also proposes a second order approximation technique which is 
interesting but requires second derivatives of the objective function and constraints. The 
second order method solves a quadratic approximating problem for a specified value of the 
parameter. The second order method will give good results in a larger region about the 
optimum than first order methods and does not appear to be as sensitive to changes in the 
active set as other methods are. However, there is still the problem of obtaining the 
Hessians of the objective function and constraints and solving the quadratic approximating 
problem. Vanderplaats and Cai (1987) feel that the second order EDS algorithm is the best 
option short of reoptimizing the problem for estimating sensitivities. But they caution that 
the method should not always be used because of its high computational cost. 

2.1.4. Variable Sensitivities 

McKeown (1980 a,c) has developed sensitivity analysis techniques for determining 
the sensitivity of design variables subject to perturbations about the optimum. This 
technique is based on an eigenvector analysis of the reduced Hessian matrix which applies 
to a variant of our standard problem (1.1)-(1.6) where no problem parameters exist . For 
unconstrained problems the major eigenvector will point in the direction of maximum 
increase of the objective function, the minor eigenvector will point in the direction of 
minimum increase of the objective function. For constrained problems the directions are 

7 


projected on the active constraints. This type of information may be useful for setting 
tolerances on design variables. 

For McKeown's algorithm, the Hessian of the Lagrangian is needed but the 
analysis is performed using only the reduced Hessian of the Lagrangian. An algorithm is 
provided for reducing the Hessian. If the Hessian is to be evaluated numerically, an 
algorithm is provided for the calculation of the reduced Hessian of the Lagrangian directly. 
This will reduce the number of extra function evaluations that are needed to conduct the 
sensitivity analysis. 

2.1.5. Other Work 

Garcia and Zangwill (1981) describe a Homotopy approach that can be used to 
solve nonlinear programming problems. They state that this approach can also be used to 
solve parametric nonlinear programming problems and is closely related to sensitivity and 
perturbation analysis. Komija and Hirabay (1984) discuss some theoretical topics involved 
in using a Homotopy approach to calculate parameter sensitivities when the active set of 
constraints changes. 

Dinkel and Kochenberger and Wong (1983) have developed an incremental 
approach for solving for the sensitivities of geometric programming problems. The 
approach is to ask the user for the new value of the parameter and then make several steps 
with corrections to reach that point. They found the smaller the step they used the more 
accurate the solution would be. 

Jittomtrum (1984) examines solving for the sensitivity of degenerate optimum 
points using the Kuhn-Tucker sensitivity equations. He provides a way to solve these 
problems using directional derivatives which provides different answers for both positive 
and negative perturbations in the parameters. Other theoretical issues for the use of 
directional derivatives to calculate optimum parameter sensitivities have been addressed by 
Janin (1984), Gauvin and Dubeau (1983), and Rockafellar, R. T. (1984). 

Zolezzi (1985) examines the conditions under which the Lagrange multipliers are 
continuous under perturbations in the problem data. This is important because Kuhn- 
Tucker sensitivity analysis uses Lagrange multipliers and rates of change of the Lagrange 
multipliers to predict the rate of change of the objective function. Comet and Laroque 
(1987) establish conditions under which the values of the Lagrange multipliers are 
Lipschitz continuous for perturbations in the problem data. 

Ganesh and Biegler (1987) have developed a sensitivity analysis based on the 
reduced Hessian. The reduction is conducted by using the equality constraints and the 
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implicit function theorem to reduce the dimensionality of the Hessian matrix that needs to 
be calculated. Their method is beneficial when there are equality constraints present in the 
formulation of the problem, because they have reduced the number of function evaluations 
required to find the required second order information numerically. Their method does not 
provide dv/9p without calculating the full Hessian of the Lagrangian. 

Rao (1987 a) and Guang-Yaun and Wen-Quan (1985) have studied the problem of 
dealing with fuzzy constraints and fuzzy objective functions. In their work they first solve 
a crisp problem then they attempt to calculate how far they can relax constraints while 
improving the objective function. To use their technique the user is required to specify 
how much violation is allowed in the constraints. Templeman (1987) reports using fuzzy 
set theory and optimization to design structures and deal with uncertainties in the problem. 

Sandgren, Gim and Ragsdell (1985) describe a problem formulation that can be 
used to obtain optimum designs with a minimum sensitivity to uncontrollable parameters. 
Their approach does not use post optimality analysis but uses a modified objective function 
to deal with the uncertainties in the problem parameters. 

The area of calculating sensitivity derivatives with respect to design variables ( i. e. 
the calculation of gradients of functions) has been an area of active research. This can lead 
to significant savings over using finite differencing. The structural optimization community 
now widely uses sensitivity analysis when the finite element method is used to analyze a 
structure. An excellent survey article of methods of sensitivity analysis for structural 
optimization is provided by Adelman and Haftka (1986). 

Haug and Arora, et al. (1977,1979,1981) have developed ways to calculate the 
gradients analytically for many structural and dynamic applications. Many of these 
methods are described in the book by Haug, Komkov and Choi (1985). 

Sobieski, et al. (1981,1982,1983,1984,1985,1986,1987) has been working on 
developing sensitivity techniques for use with multi-level decomposition techniques. 
Decomposition methods break the solution of a large problem into a system level problem 
and a group of subproblems. Each subproblem is solved using a special formulation and 
inputs from the system level problem. A sensitivity analysis is performed on the 
subproblem and the results are feed as input to the system level problem. The system level 
problem gathers all the sensitivities of the subproblems and then based on these inputs and 
others, determines the next iteration of the process. Usually, the equations (2.6) - (2.9) are 
used at the subsystem level to determine the required sensitivities, but some difficulties 
have been encountered when changes in the active set occur. 
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Schmit and Chang (1984) have developed an extension of Sobieski’s work and 
derived sensitivity equations for structural optimization problems. They derived more 
restrictive limits on the allowable perturbations than those provided by Sobieski. They 
have assumed that second derivatives of the constraints are available which is true of many 
structural problems but may not be true for other application areas. 

Schmit and Chang formulated their structural optimization problem using reciprocal 
variables and solved for the sensitivity of the dual problem. For their structural problems, 
the Hessian of the Lagrangian was diagonally dominate and the Hessian of the objective 
function was analytically available. For this class of problems good results can be expected 
even if the Hessian of the Lagrangian is inaccurate. 

Buys and Gonin (1977) developed and implemented a sensitivity analysis 
procedure for an augmented Lagrangian (AL) type code, VF01A. Their implementation is 
encouraging because they make use of the approximations of the Hessian of the Lagrangian 
that were calculated during the solution of the original problem. The results that they 
obtained using the approximate matrices were in very close agreement of those obtained by 
using the exact matrices. 

McKeown (1980 b) derives both the first and second order Kuhn-Tucker parameter 
sensitivity equations. He also provides a discussion of Fiacco's sensitivity for SUMT 
penalty functions versus Buys and Gonin's sensitivity for AL penalty functions. He 
concludes that using sensitivity for AL penalty functions should be superior to sensitivity 
by SUMT because AL produces better conditioned matrices. 

2.2. PREVIOUS WORK IN ESTIMATING PARAMETER SENSITIVITIES FOR 
CHANGES IN THE ACTIVE SET 

When the active set of constraints changes, one of the underlying assumptions 
made in deriving the Kuhn-Tucker sensitivity equations is violated. This can result in 
inaccuracies in any extrapolations based on these sensitivities since, in general, a change in 
the active constraints will result in a different set of sensitivities. Accurate sensitivity 
analysis in the presence of active set changes is also very important for efficient 
convergence of the multi-level decomposition techniques proposed by Sobieski and, in 
general, for an accurate representation of the local sensitivities. 

In the following subsections, we will first discuss the different cases that occur as a 
result of a constraint entering or leaving the active set, what effects these cases have on 
sensitivity analysis, and how changes in the active set can be predicted. We will then 
present examples of the sensitivities for the different cases which will also serve to indicate 
how the different sensitivity algorithms perform. 
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2.2.1. Cases to Consider 


When a new constraint enters the active set, or a currently active constraint leaves 
the active set, we can expect a change in the sensitivity derivatives. However, it is also 
possible that the linear independence of the constraint gradients can also be affected. For 
the discussion that follows, we define the following four cases that can result from changes 
in the active set, 

1. A constraint enters the active set and the constraint gradients are linearly 
independent. 

2. A constraint leaves the active set and the constraint gradients are linearly 
independent. 

3. A constraint enters the active set replacing an active constraint and the constraint 
gradients are linearly dependent. 

4. A constraint enters the active set and feasible region disappears. 

For Cases 1 and 2, we can expect discontinuities in the following derivatives when 
the active set changes: d2f*/dp2, 8x*/9p, and 5u*/dp. 

Case 3 is characterized by a discontinuity in the Lagrange multiplier estimates which 
causes a discontinuity in df*/dp. Since the active set changes there will also be a 
discontinuity in dx*/dp. At the point where the constraints become linearly dependent , the 
Kuhn-Tucker sensitivity equations become singular. Often what is happening for Case 3 is 
that an exchange of constraints in the active set is about to take place (i.e. the new 
constraint may replace one of the constraints that is already in the active set ). If the 
problem is not poorly formulated, we will find ourselves moving through the degenerate 
point as p increases or decreases and one of the constraints will be dropped from the active 
set. 


Case 4 is characterized as a point from which p can only be perturbed in one 
direction. If p is perturbed in the wrong direction this will cause there to be no feasible 
region and there will be no solution for the optimization problem with this value of p. Thus 
we can only perturb p in the one direction that causes the optimum path to move into the 
feasible region, and there will only exist a directional derivative for the problem in that 
direction. Case 4 can be thought of as an overconstrained design where the designer 
adjusted a parameter to the point where the design is no longer able to meet specifications. 

2.2.2. Prediction of when the Active set will Change 

Barthelemy and Sobieski (1983 a) have observed that the accuracy of extrapolations 
of the objective function deteriorates rapidly when the active set changes. From section 

2.1.2, we saw that we can use equations 2.10 and 2.1 1 to predict where the active set will 
change, thus we can use this information to predict when the extrapolations will deteriorate. 
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A problem with bounding Ap by equations 2.10 and 2.1 1 is that the estimate is only 
good for the first constraint that is encountered because once the active set changes the 
search direction to the new optimum will change (the discontinuity in dx/dp). Thus, it 
becomes very difficult to estimate when or which constraint will leave/enter the active set 
second. This problem will be addressed in section 6. 

The merit of using equations 2.10 and 2.1 1 to predict when the active set will 
change was discussed by Adelman and Haftka (1986). They state, "The effectiveness of 
using this approach (equations 2.10 and 2.1 1) is still in doubt with positive results being 
obtained by Schmit and Chang (1984) and negative results being obtained by Barthelemy 
and Sobieski (1983 a)". We feel that the positive results that were obtained by Schmit and 
Chang are due to problem linearity and the changes in the active set that they encountered 
being case 1 and case 2 changes. We feel that the negative results obtained by Barthelemy 
and Sobieski are due to nonlinearity of the problem and also a case 3 change in the active 
set taking place. As we will see later in this report, the consequences on sensitivity 
derivatives of case 3 changes in the active set are often much more severe than case 1 and 
case 2 changes. 

2.2.3. An Example of Case 1 and 2 


The effect of a constraint entering or leaving the active set (Cases 1 and 2) can best 
be demonstrated by a simple example from Vanderplaats and Yoshida (1985). 

Minimize f(x) = 2xi2 - 2xip + p2 + 4xi - 4p 
subject to: gi = 4p + xi > 0 

The Lagrangian will be 

L(x,u) = 2xi2 - 2xip + p2 + 4xi -4p - ui(4p + xj) 
for p = 0, the optimum is f(x*) = 0, xi* = 0, gi = 0, and ui = 4. 

This example will illustrate a constraint leaving the active set (case 2) as p increases. 
The same example can be used to illustrate a constraint entering the active set (case 1) if we 
use a different starting value of p. 

To demonstrate the methods we have talked about, we will calculate the sensitivity 
estimates using four representative methods: the first and second order Kuhn-Tucker 
method, and first and second order extended design space method. We will conclude with 
a comparison of the various methods used to solve the problem. 

To solve for the sensitivity by Kuhn-Tucker equations we use equation (2.6) to 
provide the following system of equations 


(2.13) 

(2.14) 

(2.15) 
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which yield 


dxf 

L'SiTJ 


■ [4] 


(2.16) 


|r =- 4 < 2 - 17 > 

2^.= -18 (2.18) 

From equation (2.7) we can determine the sensitivity of the objective function with respect 
to the parameter p to be, 

3F = |r + Jrr T =-4 + 4 (-4) = -20 (2.19) 

The active set will change when the Lagrange multiplier of the constraint goes to 
zero, which can be estimated by equation (2.10) 


Ap = 


- ui 


-4 


0UA -18 

W) 


= 0.2222 


( 2 . 20 ) 


therefore we are assured of reasonable results for extrapolations for which Ap less than 

0 . 2222 . 


For example, a linear approximation by equation (2.8) to estimate the value of the 
new optimum produce 

df f\ . A / /, A\ AA 1 (221) 


f new = f*+ Ap^ = 0 + Ap(-20) = -20Ap 


A quadratic estimate of the new value of the objective function can be made by 
evaluating the following equation found in Fiacco (1983), McKeown (1980 b), and 
Sobieski and Barthelemy (1983) 

d 2 f 0 2 L d 2 L 3x1 Oui 9gi n 

dp 2 3p 2 dxjdp dp dp op ' ’ 

v/hich produces d 2 f/dp 2 = 82. Using the quadratic estimate for the value of the objective 

function we obtain 


fnew = f* + Ap^ + 0.5Ap||Ap = -20Ap + 41 Ap2 


(2.23) 


The same predictions can be made by Vanderplaats' extended design space 
algorithm. We begin by formulating the following direction finding problem for decreasing 
values of p, where X2 represents the parameter p, and X3 is an additional variable to ensure 
that p has the required sign. 


minimize 4xi - 4x2 - c X3 


(2.24) 


13 


subject to: xi + 4 x 2 ^ 0 

-X2 - X3 > 0 
1- (xi2 + x22)>0 


(2.25) 

(2.26) 
(2.27) 


For c = 1000, the solution is xi = .970142, X 2 = -.242536, X 3 = .242536 
yields the following estimates of the sensitivity derivatives 


£=-2° 

dp 

d*L__ 4 
dp - 4 


which 

(2.28) 

(2.29) 


For increasing values of p we obtain the following subproblem 


minimize 4xi - 4 x 2 - c X 3 (2.30) 

subject to: x i + 4x2>0 (2.31) 

x 2 - x 3 > 0 (2.32) 

1- (xi 2 + X 2 2 )^ 0 (2.33) 


When this problem is solved, the resulting sensitivities are sensitive to the value of 
the parameter c. The solution for several values of c are presented below in Table 2.1. 


Table 2.1 The effect of "c" on EDS sensitivity 


variable 

c =1000 

c=500 

c =100 

c =10 

c= 1.0 

O 

II 

O 

b 

xi 

-0.398E-2 

-0.79E-2 

-0.388E-1 

-0.2763 

-0.624 

-0.707 

x 2 

0.99999 

0.99996 

0.99924 

0.9611 

0.9611 

0.707 

x 3 

0.99999 

0.99996 

0.99924 

0.9611 

0.9611 

0.707 

df/dp 

-4.016 

-4.0316 

-4.155 

-5.150 

-7.196 

- 8.0 


From this table it is clear that the choice of c will effect the sensitivity derivatives. For 
demonstration purposes c = 10 was chosen, this yielded the following sensitivity 
derivatives. 

§=•5.1502 (2.34) 

|j|-= -.28756 ( 2 . 35 ) 

Vanderplaats and Yoshida (1985) report that the value of c has little effect on the 
EDS algorithm. However Vanderplaats and Cai (1987) report that after further research the 
value of c will effect the accuracy of the EDS procedure. 

Using Vanderplaats second order extended design space algorithm provides exact 
answers for the sensitivity for this problem. 
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Figure 2.1 illustrates the accuracy of various methods. We can see that when the 
active set changes at Ap = 0.222 the predictions become less accurate. 



g(l) leaves p 

the active set 


Figure 2.1 A plot of various optimal values of f with respect to p 

Figure 2.2 illustrates the location of the optimum value of xi as a function of p, as 
predicted by various algorithms. When the active set changes there is a discontinuity in the 
rate of change of the optimum value of xi with respect to p (i. e. dxi/dp is discontinuous at 
the point). 




I 


Figure 2.2 A plot of various optimal values of xj with respect to p 


From figures 2. 1 and 2.2 it is possible to draw some conclusions about the relative 
performance of the four different methods that were used to obtain sensitivity information. 
Using the first order Kuhn - Tucker method we see that the solution follows the inequality 
constraint in both the positive and negative direction. The linear estimate of the new value 
of xj is accurate for small changes in p less than 0.2222. But for values of p greater than 
0.2222, the active set has changed and large errors in the predictions are introduced. This 
is also true for the linear prediction for the value of the objective function. 

The second order Kuhn - Tucker estimate of the value of the objective function is in 
exact agreement in the region where the active set remains the same, as seen in figure 2.1. 
However after the active set changes the predicted value of the objective function is a poor 
predictor of the actual value of the optimum. 

The first order extended design space provides the same results as the first order 
Kuhn-Tucker sensitivity for decreasing values of p. For increasing values of p we see that 
the search direction changes. This approximation appears to overcome the constraint 
leaving the active set, but it is a poor predictor of the actual value of the optimum for small 
variations in p. For other values of the parameter "c" we will obtain similar values for the 
sensitivity derivatives. 

The second order extended design space provides the exact values of the locations 
of the optimum value of the objective function. This is because the approximating problem 
that is formulated is the same as the original problem. 

With this simple example we have demonstrated the effect of a constraint leaving 
the active set on the algorithms for estimating parameter sensitivity. We can see from this 
example that, as we might anticipate, using second order estimates can produce more 
accurate extrapolations. In fact, only the second order extended design space algorithm 
provided good results after the constraint left the active set. However its usefulness is 
diminished by the need for second derivatives which can be computationally expensive to 
obtain. 

2.2,4. Example of Case 3 

Recall, that Case 3 is characterized by the adding of a new constraint to the active 

set and the gradients of the active constraints become linearly dependent. When the 

gradients of the constraints are linearly dependent the Lagrange multipliers will not be 

uniquely determined and the Kuhn-Tucker optimality conditions cannot be uniquely 
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verified. This also results in a discontinuity in the Lagrange multiplier sensitivities. 

When the constraint gradients become linearly dependent for a value of one of the 
parameters it is assumed that this is only a temporary condition. If the user is interested in 
the effect of changing the parameter on the optimum then this information can be obtained 
on either side of the singular point. 


This behavior is demonstrated in the following example 

minimize: f = xj2 + (P - 1)2 (2.36) 

subject to: gi = 3 xi + 2 P - 10 > 0 (2.37) 

g2 = 2 xi + 3 P - 10 > 0 (2.38) 

When P = 2, the minimum f* = 5 occurs at xi* = 2. At this point , both constraints are 
active, and the gradients of the constraints are not linearly independent. The Lagrange 
multipliers will be in the family 

ui,u2€ {3 ui + 2 U2 = 4, ui > 0, U2 > 0} (3.39) 


At this point, df*/dp, 9x*/3p and 3u/3p can not be uniquely determined. Results 
for these derivatives can be developed if we consider positive and negative changes in p 
separately on either side of this degenerate point which we shall indicate by 3x/3p+ for 
increasing values of p and 3x/3p- for decreasing values of p. 

Figure 2.3 presents the sensitivity plots for this problem. Figure 2.3 (a) and (b) 
represent the first order predictions of the new values of the Lagrange multipliers for this 
problem. For this problem the linear predictions agree with the optimum Lagrange 
multipliers. There is a discontinuity at Ap = 0.0, therefore there will only be directional 

derivatives for these values. Figure 2.3 (c) represents linear predictions of the new value 
of the objective function. Notice again that there is a discontinuity in the slope of the 
prediction and we can not determine df* /dp for Ap = 0. Therefore df*/dp will not exist for 
this value of p. Figure 2.3 (d) represents the predicted location of xi and we notice the 
same situation as we have for df* /dp. 
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Figure 2.3. A comparison of the Sensitivity of a problem with a Linear Dependence in the 

Constraint normals 

2.3. SUMMARY 

Sensitivity analysis is now routinely used in linear programming (Falk and Fiacco 
1982) and most linear programming algorithms provide modules for the calculation of 
sensitivities. This has not been the case for applications of nonlinear programming. The 
most common use of sensitivity derivatives has been in the area of structural optimization 
and in work done for decomposition methods. Some of the reasons for this may be due to 
a lack of understanding about how to perform sensitivity analysis for nonlinear problems, 
or to a lack of established procedures and supporting software that make the analysis more 
readily available to the average user. The largest contributor to its lack of use is probably 
the difficulty involved in implementing the current theory and methods. 

An assessment of the methods discussed in Section 2. 1 and demonstrated in the 
examples in Section 2.2.4 leads to the following conclusions about the current state of the 
art of parameter sensitivity analysis: 

1 . The Kuhn-Tucker sensitivity equations (2.6) accurately define the desired 
sensitivities assuming no changes in the active constraints. To implement these 
equations, however, requires second order information about the Hessian of the 
Lagrangian, and the change in the gradient of the Lagrangian with respect to the 
parameter. Both of which are difficult to obtain reliably for all but a few special 
cases. 


18 



2. The Extended Design Space (EDS) method provides sensitivity information 
without the need for the second order information required of the Kuhn-Tucker 
method. However, the sensitivity estimates are effected by a choice of a 
formulating parameter c, and may not give the same directions as those obtained 
from the Kuhn-Tucker method. 

3 . Changes in the active constraint set will effect the accuracy of any of the 
methods and may limit the region upon which extrapolations to the design can 
be relied upon. 
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3. 


New method for Estimating Parameter Sensitivity 


In this chapter a new method for estimating parameter sensitivities based on the 
Recursive Quadratic Programming(RQP) method is described. We begin with a brief 
description of the RQP method and the advantages it provides for estimating sensitivities. 
Next, we present the RQP based algorithm for estimating parameter sensitivity that exploits 
the advantages of the RQP method discussed in the previous section. This is followed by a 
comparison of the new method with existing methods based on the type of information that 
is being produced and the number of function evaluations required. Finally, a discussion 
is presented of potential problems that may be encountered with the new RQP sensitivity 
method. 

3.1. RQP METHODS 

The RQP method has been on the forefront of recent research in optimization 
algorithms and has been emerging as one of the most efficient methods available for 
solving small to medium sized, general nonlinear programming problems (equations 1.1- 
1.6). State of the art RQP methods have been developed by many researchers, such as, 
Powell (1983), Schittkowski (1984), Gill, Murray and Wright (1986) and Bartholomew- 
Biggs (1986,1987) to name a few. The algorithm has been tested against other general 
nonlinear programming algorithms by Schittkowski (1980), Ecker and Kupferschmid 
(1984), Belegundu and Arora (1985). The results of these tests have shown the RQP 
method to be one of the most efficient algorithms available for the solution of nonlinear 
programming problems. 

All RQP methods use the same basic strategy of linearizing the constraints and 
approximating the Hessian of the Lagrangian to form a quadratic programming (QP) 
subproblem. The QP subproblem is then solved for the search direction, s, and a new 
estimate the Lagrange multipliers of the constraints. The QP subproblem has the form 

Minimize l/2sTBs + s T Vf (3.1) 

subject to VhTs + h = 0 (3.2) 

VgT s + g >0 (3.3) 


where B is an approximation to the Hessian of the Lagrangian which is normally 
constructed by variable metric methods. The Lagrange multipliers of the constraints for the 
original problem (equations 1.1 -1.6) are estimated by the Lagrange multipliers of the 
constraints in the QP subproblem (equations 3. 1-3.3). The search direction s is then used 
to calculate a new estimate of the optimum 
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xit+i = xit + as (3.4) 

where a is determined by minimizing a line search penalty function P of the following 
general form, 

P(x,u,v,R) = f(x) + R*Q(h,g,u,v) (3.5) 

where Q represents some combination of the constraints and the Lagrange multipliers. The 
penalty function attempts to assure that both the objective function and the violation of the 
constraints are reduced. As the method converges, the optimal step length a generally 
approaches 1. 

RQOPT, a typical RQP algorithm, was used in our research. A summary of the 
algorithm that is used by RQOPT is presented here, a full description of RQOPT can be 
found in the users manual (Beltracchi and Gabriele 1987 a), or Beltracchi (1985), Gabriele 
and Beltracchi (1986,1987 b). There were several modifications that were made to RQOPT 
for this work and these will be discussed in section 4.1 of this report. 



Figure 3.1 Flow Chart for RQOPT 

Figure 3.1 shows the basic steps that are used by the RQOPT program. The 
RQOPT algorithm begins with an initial estimate of the location of the optimum and several 
algorithm parameters that have been set by the user. The first step of the algorithm is to 
identify the active constraints, it is important that the proper constraints are chosen to be in 

the active set as this can effect the rate of convergence of the algorithm and, for our 
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purposes, the approximation of the Hessian of the Lagrangian. Algorithm parameters are 
available to allow the user to control which constraints are considered active during the 
course of the optimization. 

The second step is to calculate the gradients of the objective function and the 
constraints that are in the active set and then update the approximation of the Hessian of the 
Lagrangian. The update of the Hessian is performed using the BFS variable metric update 
with modifications specified by Powell (1977). 

The third step is to solve a quadratic programming subproblem (equations 3. 1-3.3). 
The QP subproblems generated by RQOPT are solved by OPTQP, a special implementation 
of the reduced gradient method. If the subproblem has no feasible solution, the active set is 
redefined by dropping constraints from the active set until a feasible subproblem can be 
found. 


The line search for the next point x it+1 makes up the fourth and fifth steps of the 
algorithm. An initial step size for the line search is determined in the fourth step such that 
constraints not in the active set are not excessively violated. The line search is performed in 
the fifth step, and if a step of a = 1 satisfies the line search criteria, then that step is taken 
and the line search ended. 

The sixth step updates the penalty parameters used in the line search, and the 
Lagrange multiplier estimates. We then return to start another iteration. 

There have been several different variants of the RQP method proposed. Some of 
the variants are discussed in Beltracchi (1985). The major differences in RQP algorithms 
are in the form of the line search objective function (equation 3.5) and the formulations of 
the QP subproblem (equations 3. 1-3.3) that are used. Research continues on these areas 
but no one formulation has yet to prove itself clearly superior. 

Some of the penalty functions that have been proposed for (3.5) are a h exact 
penalty function (Fletcher 1984, Powell 1987), a I 2 quadratic loss penalty function 
(Bartholomew-Biggs 1980) or an augmented Lagrangian (Chen, Kong and Cha 
1987,Bartholomew-Biggs 1985, 1987). The penalty function's parameters are adjusted 
after each iteration, and how the parameters are updated effects the convergence of the 
method. 

There are two basic philosophies for forming the QP subproblem for RQP 

methods, the inequality constrained (IQP) formulation and equality constrained (EQP) 

formulation. The most common is the IQP approach which uses a subproblem of the form 

of equations 3. 1-3.3. The EQP approach linearizes only a subset of the inequality 
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constraints and considers these as equality constraints in the subproblem(i.e. equation 3.3 
is considered to be an equality constraint). A discussion of the advantages and 
disadvantages of the IQP and EQP subproblem formulation can be found in (Bartholomew- 
Biggs 1987,1986,1982, Zhou and Mayne 1985, Schittkowski 1983, Murray and Wright 
1982, or Powell 1978). 

Although the method does perform well, it does have some disadvantages. In 
general, the method produces a series of infeasible points while approaching the solution 
which may pose a problem for some problem formulations. RQP methods are also 
sensitive to variable and objective function scaling and no good scaling algorithms have 
been proposed. Finally, the best penalty function or algorithm for updating the penalty 
parameters for the line search is still a subject of a great deal of research in these methods. 

On the plus side, the following advantages have been attributed to the method. In 
terms of number of function evaluations, this method appears to be one of the most 
efficient methods available. This has been demonstrated in any of the published 
comparison studies in which codes for these methods participants. The method does not 
require a feasible starting point which means there is no special phase 1 search employed as 
in the GRG method or the feasible direction method. Although, as mentioned above, the 
method is sensitive to variable and objective function scaling, it is not sensitive to 
constraint scaling. Finally, the RQP method provides an estimate of the Hessian of the 
Lagrangian, which can be useful for other purposes, and it is very efficient at locating an 
optimum when the starting point is close to the true optimum. Both of these last 
advantages will be exploited in the next section which describes a method for sensitivity 
estimation based on the RQP method. 

3.2. PROPOSED ALGORITHM FOR PARAMETER SENSITIVITY 

In reviewing the current methods for sensitivity analysis in chapter 2, we recall that 
to employ the Kuhn-Tucker sensitivity equations required second order information about 
the Lagrangian. For most engineering problems this type of information is often not 
available in closed form, and estimation techniques would be prone to truncation and 
numerical errors. Therefore, the application of these equations to a broad spectrum of 
engineering applications is limited. 

One proposal mentioned in chapter 2 to circumvent these problems was suggested 
by Armacost and Fiacco (1974) and McKeown (1980 b). Their proposal to estimate the 
sensitivities without estimating the higher order information was given in equations 2.4 and 
2.5, 
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df* _ f(x*,p + Ap) - f(x*, p - Ap) 
d P “ 2 Ap 

9x* _ x*(p + Ap) - x*(p - Ap) 
d P “ 2 Ap 

These equations represent the use of differencing techniques to estimate the sensitivities, 
where the values f(x*,p + Ap), x*(p + Ap), etc. are determined by reoptimizing the 
problem for the new values of the parameter. For most algorithms, particularly penalty 
function based methods, the reoptimizations would be a non-trivial task requiring a 
considerable number of function evaluations. However, this is the type of problem where 
the RQP method is considered to be very effective. The goal of the new algorithm is to 
exploit the strengths of the RQP method to estimate sensitivities by these differencing 
techniques. 

The RQP method possesses two characteristics that we felt can be exploited for 
determining parameter sensitivities: (1) an approximation to the Hessian of the Lagrangian 
is developed, and (2) if this approximation is exact (or close) then the RQP method quickly 
and efficiently solve the reoptimization problem used in the difference equations. 
Essentially, if we can develop good Hessian approximations, the RQP method is equivalent 
to applying Newton's method to solve the Kuhn-Tucker conditions for the perturbed 
problems which should require only 1 or 2 iterations of RQP 1 . The small number of 
iterations, coupled with the fact that the RQP method should require only a one step line 
search, should allow the reoptimizations to occur without the need for many function 
evaluations. 

Based on the above arguments, we propose the following procedure to calculate 
parameter sensitivity derivatives (for cases where there are no changes in the active set for 
small variations in the paramters 2 ). 

Step 0. Given an optimal solution x*, f*, u*, an active set of constraints, and an 
approximation to the Hessian of the Lagrangian, all achieved by convergence 
of the RQP method. 

(the * notation is used to denote optimum values) 

Step 1. Perturb the fixed parameter pi to pi+ = pi 0 + Apj where Apj is some small 
perturbation to pi 

Step 2. Perform one complete iteration of the RQP method to find: 
f + the estimated value of the optimum objective function 
x + the estimated value of the optimal of the design variables 

1 We can expect only one or two iterations of RQP if we can adequately approximate the perturbed problem 
with a quadratic function. Due to the small region of interest, a quadratic approximation should be good. 

2 At points where the active set changes then modifications discussed in chapter 6 must be used to calculate 
directional derivatives. 
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u + the estimated value of the optimum Lagrange multipliers 
gj + j £ Active set 

(as predicted by the RQP method for pi = pi+) 

Step 3. Perturb the fixed parameter p, to pf = pj° - Apj 

Step 4. Perform one complete iteration of the RQP method to find: 
f" the estimated value of the optimum objective function 
x- the estimated value of the optimal of the design variables 
u' the estimated value of the optimum Lagrange multipliers 
gj~ j € Active set 

(as predicted by the RQP method for pi = pf) 

Step 5. Obtain estimates for the sensitivity derivatives from the following central 
difference approximations 


df* f+ - f- 

(3.6) 

d P 2Ap 

dx* x + - x" 

(3.7) 

"^P - " 2Ap 

3u* u + - u* 

(3.8) 

cx, 

<1 

<N 

II 

k 


Step 6. Estimate the sensitivity of the inactive constraints by 

= JLl — lJl_ j £ Active set (3.9) 

d P 2Ap J 


In addition to the algorithm described above, the following variants of the basic 
algorithm are also proposed 


1 . Forward differencing, For this variant we would omit steps 3 and 4 and then 
use a forward difference approximation (equation 3.10 instead of equations 3.6- 
3.9) to approximate the derivatives 

=- 3 - + ~ a (3.10) 

d P Ap 

where q can represent f*, x, u, and the inactive constraints. We may want to 
use this formulation because it requires less function evaluations than the central 
difference approximation. However, the forward difference approximation is 
more susceptible to roundoff and truncation errors and requires a more accurate 
optimum to yield good sensitivity derivatives. 

2. Forward differencing using 2 iterations of the RQP method. This variant is 
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similar to variant 1, but we would perform 2 iterations of RQOPT in step 2. This 
will yield a more accurate estimate of the optimum of the perturbed problem. 
When we use this option we can also update the approximation to the Hessian of 
the Lagrangian, or adjust the perturbation Ap to obtain a more accurate estimate 
of the derivatives. 

3. Central differencing using 2 iterations of the RQP method. This variant would 
perform two iterations of RQOPT in steps 2 and 4 of the basic algorithm. As in 
variant 2 we can update the Hessian approximation during each iteration or adjust 
the perturbation Ap to obtain a more accurate estimate of the derivatives. This 
variant is the most computationally expensive of the proposed variants. 

When there are many parameters that the user needs to obtain sensitivities for then 
the user may want to use variant 2 or variant 3 to calculate the sensitivities for the first few 
parameters. This will allow a more accurate estimate of the Hessian of the Lagrangian to be 
constructed. After an accurate estimate of the Hessian of the Lagrangian is built, the user 
should switch to either the baseline or variant 1 to obtain the sensitivities of the remaining 
parameters. The Kuhn-Tucker sensitivity equations may also be used with the Hessian 
approximation, after a good estimate of the Hessian of the Lagrangian is built. However 
the Kuhn-Tucker sensitivity equations also require 3V x L/dp be calculated and this term 
may be subject to numerical noise because V X L = 0. 

3.3. COMPARISON TO EXISTING METHODS 

This section provides a derivation that indicates the performance that is expected 
from the new sensitivity algorithm. This section also presents a comparison between the 
RQP based method and two existing methods described in chapter 2 based on the number 
of function evaluations required to estimate the sensitivities. 

3.3.1 . Demonstration of Equivalence of New Method to Kuhn-Tucker Method 

This section will show that the finite difference approximations obtained by the 
proposed method are in fact equivalent to the sensitivities obtained by solving a modified 
set of Kuhn-Tucker sensitivity equations. The modification of the Kuhn-Tucker sensitivity 
equations involves replacing the Hessian of the Lagrangian with the approximation B, 
obtained from the RQP method. 

The following assumptions are made for this derivation; no equality constraints are 
present, the base optimal point is stable 3 , and the gradients are continuous. The derivation 

3 A stable point is defined as a point where the acitve set does not change for small variations in the 
parameters 
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in the presence of equality constraints does not change too much but the equality constraints 
were left out to simplify the notation. If the base point is not stable then this derivation can 
be used to find directional derivatives; this will be discussed at the end of this section. If 
the gradients are not continuous then we cannot even be assured of an optimum point since 
the assumption of continuity is also made for the derivation by the Kuhn-Tucker method. 


We begin by restating the Kuhn-Tucker Sensitivity equations 


' V^L -V x g~ 

5x~ 

¥ 

3u 

+ 

“av x L _ 

"3T 

V x g T 0 . 

L5pJ 
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We strive in this derivation to show that the proposed method is equivalent to estimating the 
sensitivities using modified version of equation (3.1 1) that replaces V^L with B obtained 

from the RQP method. If this is the case, then we can anticipate the kind of accuracy to 
expect and where the possible sources of error will result. 


If we examine the equations (3.6-3. 8), used by the proposed RQP sensitivity 
method we see that these provide finite difference approximations to the sensitivity 
derivatives of the objective function, design variables, and Lagrange multipliers with 
respect to pj. The derivatives are defined by the following 

df* = lim | T*(x*+Ax,pO+Ap) - f*(x*,pO) 
dp Ap-»0 y Ap 


(3.12) 


9x* lim fx*(p°+Ap) - x*(p°)A 

■5T = Ap-»0 Ap J 


3x* 

Ax=^-Ap 


8u* lim / ii*(p()+Ap) - u*(pO)> 
~3p” - A P-*° Ap , 

where p° represents our base point. 


(3.13) 

(3.14) 

(3.15) 


The RQP subproblem for the simplified case where the active constraints remain 
active and there are no equality constraints can be written as 

min 1/2 s T Bs + s T V x f (3.16) 

subject to s T V x gj + gj = 0 j e Active Set (3.17) 

where B is the approximation to the Hessian of the Lagrangian and the inequality 
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constraints gj are considered as equality constraints. 


If we assume 4 a step length of a = 1 is used in the line search (equation 3.4) we 
can rewrite equation 3.4 in terms of x' - x as 

s = x' - x (3.18) 

where x' is the new estimate of x*. Substituting equation 3.18 into equations 3.16 and 
3.17 we obtain the following subproblem which is minimized with (x* - x) as the design 
variables 

min 1/2 (x’-x) T B(x'-x) + (x’-x) T V x f(x,p + ApO (3.19) 

subject to (x'-x) T V x gj(x,p + ApO + gj(x,p + ApO = 0 j e Active Set (3.20) 
We can now state the optimality conditions for the subproblem represented in equations 
(3.19-20) as 

B(x'-x) + V x f(x,p + ApO - u'V x gj(x,p + ApO = 0 (3.21) 

(x'-x) T V x gj(x,p + Apt) + gj(x,p + Apj) = 0 j e Active Set (3.22) 

Here u' represents the estimated value of the Lagrange multipliers at the new optimum. 

Now we substitute into equation (3.21) the following definitions of zero 

V x L(x,pO) = 0 = V x f(x,p0) - uV x g(x,pO) (3.23) 

uV x g(x,p°+ Ap) - uV x g(x,p°+ Ap) = 0 (3.24) 

This will yield 

B(x'-x) + V x f(x,p° + Ap) - u'V x gj(x,p°+ Ap) - (V x f(x,p°) - uV x g(x,p°))+ 
uV x g(x,p°+ Ap) - uV x g(x,p°+ Ap) = 0 (3.25) 

Rearranging we obtain 

B(x'-x) - u’V x gj(x,p°+ Ap) + uV x g(x,p°+ Ap) + 

(V x f(x,p° + Ap) - uV x g(x,p° + Ap)) - (V x f(x,p°) - uV x g(x,p°)) = 0 (3.26) 
Rearranging further and writing in terms of the Lagrangian function we obtain 

B(x'-x) - (u' - u)V x gj(x,p°+ Ap) + V x L(x,u,p°+ Ap) - V x L(x,u,p°) = 0 

(3.27) 

Now we will divide equation (3.27) by Ap and take the limit as Ap goes to zero to 
4 A common assumption for RQP methods 
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obtain 


lim B(x'-x) - (u' - u)V x gj(x,p<4 Ap) + V x L(x,u,p°+ Ap) - V x L(x,u,p°) = 
Ap— >0 Ap 

(3.28) 

Using the additive and multiplicative properties of the Limit function we obtain 


lim f 'V x L(x,u,p°+ Ap) - V x L(x,u,p°)^ 

Ap J 


(3.29) 


Now we can use the definition of a derivative of some function h with respect to some variable p 

(3.30) 

Applying the definition of 5x/9p, 9u/3p to (3.29) we obtain 

(3.31) 


9h_ lim h(p+Ap)-h(p) 

<5p Ap— »o ^ 


„dx 9u lim „ , n A N dV x L(x,u,p°) 

®3p ' 3p ip ^o v ^i <x ’P 0+ + ^- E -=° 


W 


If we use the standard assumption that the functions are twice continuously differentiable 
we can state 


a p — >o V xgj( x »P° + Ap) = V x gj(x,p°) (3.32) 

<• md now substituting equation (3.32) into equation (3.31) we obtain 

du . 3 V x L(x,u,p°) „ 

B Sf - Vjgjfx.pO)^ + gj =0 (3.33) 

The equation above (equation 3.33) represents the first part of the Kuhn-Tucker sensitivity 
equations with the approximation B instead of the Hessian of the Lagrangian. 


The next step in this derivation is to examine equation (3.22) in terms of p° + Ap 
we obtain 


(x’-x) T V x gj(x,p° + Ap) + gj(x,p° + Ap) = 0 j e Active Set (3.34) 

Now we can subtract gj(x,p) = 0 from equation (3.34) to obtain 

(x’-x) T V x gj(x,p° + Ap) + gj(x,p° + Ap) - gj(x,p°)= 0 j e Active Set (3.35) 
If we divide equation (3.35) by Ap and take the limit as Ap goes to zero we can write 
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lim ^(x , -x) T V x gj(x,p° + Ap) 
Ap->0 ^ Ap 



= 0 (3.36) 


Using the additive and multiplicative properties of the limit function we obtain 

““ V xgi (xp° + Ap)* lim ('8j ( ’‘.P° + A P> - gj <x 'P 0 ) 

Ap ^0 xgj( ,P P) Ap^ Ap J Ap-»0^ Ap 


Again using the definition of a partial derivative of (equation 3.30) and we obtain from 
equation (3.37) 

V ^j( x ’P° + A P)*^ + ^p i=0 < 3 - 38 

Using the results in equation (3.32) we obtain 

V x gj ( x , pO )*|^+^=0 0.39 

Which represents the second part of the Kuhn-Tucker sensitivity equations. 

Now equations (3.33) and (3.39) can be assembled into matrix form to yield 





(3.40) 


Equation 3.40 is the same as equation 3.19 with the exception that equation 3.40 uses the 
approximation, B, of the Hessian of the Lagrangian in place of, V^L, the true Hessian of 


the Lagrangian. Referring to (3.40) as the modified Kuhn-Tucker equations, we see that 
the proposed method is principally a difference approximation to the modified Kuhn- 
Tucker equations. This implies that if B is a good approximation of the Hessian of the 
Lagrangian, and a proper choice can be made for the difference parameter that minimizes 
truncation and roundoff errors, then we can produce sensitivity derivatives without the 
need to obtain or estimate the second derivatives required of the Kuhn-Tucker method. 


Several examples were tested to see if the sensitivity derivatives estimated by the 
RQP method with one iteration converged to the value of sensitivity derivatives estimated 
by the Kuhn-Tucker sensitivity with the approximate Hessian. From these examples we 
observed that the sensitivity derivatives delivered by the new RQP algorithm are close to 
the derivatives approximated by the Kuhn-Tucker method with the Hessian approximation. 
One of these examples is presented here to show this agreement. 
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Test problem 2 (which is described in the appendix) is used to demonstrate the 
equivalence of the new method to the Kuhn-Tucker method. The starting point of x° = 
(1.1, 1.1, 1.1) was used. RQOPT(with the BFS update and H° = I) solved the problem in 
one iteration and yielded the following approximation to the Hessian matrix 5 


Happrox — 


r3 2 2 ] 
2 3 2 
L2 2 3 - 


If we use this Hessian approximation to solve for the sensitivity of parameter 1 by equation 
(3.40) we will obtain the following system of equations 


I¥1 


r3 2 2 - 1 i 


5x2 


r-12-i 

2 3 2 -1 




5 

2 2 3 -1 


5x3 

+ 

0 

Ll 1 1 0 J 




1 . 


L|rJ 


the solution of these equations yields 

= (9.33333,-7.66666,-2.66666) 

I- 

(-4.66666) 


(3.41) 


(3.42) 

(3.43) 


The RQP based sensitivity algorithm calculated the following sensitivity derivative 
approximations. 

=(9.33333,-7.66666,-2.66666) (3.44) 

= (-4.66666) (3.45) 

The above derivatives were calculated using the RQSEN program (described in section 4 
and the appendix of this report) with a perturbation of Ap = 0.G001 (using central 
differencing, equations 3.7, 3.8) and one iteration of RQP to solve the perturbed problems. 

If the base point, p°, is unstable (degenerate) we can use a similar derivation to 
calculate directional derivatives, which will be useful for predicting the sensitivities of the 
design variables and Lagrange multipliers. The use of directional derivatives will be 
discussed in section 6. 

5 The Hessian approximation for problem 2 is not close to the true Hessian of the Lagrangian ( given in the 
appendix of this report). This is because the starting point was chosen to produce a poor approximation so 
we could clearly indicate the performance of the RQP sensitivity method in comparison to the Kuhn-Tucker 
sensitivity method with the approximate Hessian from RQOPT 
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3.3.2. Performance Comparison with Other Methods 


This section compares the RQP based method to two of the methods discussed in 
chapter 2; the Kuhn-Tucker method, and the extended design space (EDS) method. The 
comparison is based on Table 3.1 which examines the number of function evaluations 
required by each method to calculate parameter sensitivities df*/dp, 3x*/3p and du*/dp 
(assuming that when the optimum is found that the Kuhn-Tucker conditions have been 
checked, this means that V X L and the Lagrange multipliers are known before the 
sensitivity analysis is performed). It is assumed that the objective function and constraints 
are interrelated 6 . It is also assumed that problem linearity or problem form are not 
exploited in calculating parameter sensitivities. 

The first row of Table 3.1 represents the methods used in this comparison. The 
second row represents the number of function evaluations required to calculate the 
sensitivity derivatives for the first parameter. Subsequent parameters may require fewer 
evaluations for some methods. 

The first column of Table 3.1 represents the number of variables present in the 
problem. The second column represents the amount of work required to solve for the 
sensitivities using the Kuhn-Tucker sensitivity equations. The third and fourth columns 
represent the number of function evaluations required by the EDS algorithm. Column 3 
represents the first order method and column 4 the second order method. The fourth 
column, RQP 1, indicates that forward difference approximations were used to calculate the 
gradients. The fifth column RQP 2, also uses forward difference approximations but 2 
iterations of RQOPT are allowed during the reoptimization. The fifth column RQP 3 
represents the amount of work required for the base line algorithm using central difference 
approximations. The sixth column RQP 4 represents using central difference 
approximations with 2 iterations of RQOPT. 

If the objective function sensitivity is calculated by equation 2.12 
(df*/dp = 9f/dp - u dg/dp) then assuming that objective and constraint information can 
both be obtain in one call, only one extra function evaluation is required to determine df/9p 
and dg/dp. However, if one wants the design variable and Lagrange multiplier sensitivity 
then some other equations must be used. 


6 The value of the objective function and all of the constraints are calculated by one subroutine 
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n 

Kuhn-Tucker EDS (1 st ) EDS (2 nd ) 

RQP 1 

RQP 2 

RQP 3 

RQP 4 


3n 

-+-+ 1 

(n+l)2 3(2+11 
2 + 2 

n+1 

2n+2 

2n+l 

4n+2 

1 

4 

1 5 

2 

4 

3 

6 

2 

6 

1 9 

3 

6 

5 

10 

3 

10 

1 14 

4 

8 

7 

14 

4 

15 

1 20 

5 

10 

9 

18 

5 

21 

1 27 

6 

12 

11 

22 

10 

66 

1 77 

11 

22 

21 

42 

15 

136 

1 152 

16 

32 

31 

62 

20 

231 

1 252 

21 

42 

41 

82 

40 

861 

1 902 

41 

82 

81 

162 


Note: RQP 1 uses forward difference approximations and one iteration to solve the perturbed problem 

RQP 2 uses forward difference approximations and two iterations to solve the perturbed problem 
RQP 3 uses central difference approximations and one iteration to solve the perturbed problem 
RQP 4 uses central difference approximations and two iterations to solve the perturbed problem 


Table 3.1 Comparison of Various Algorithms for use in sensitivity analysis 


The following observations can be drawn from this table. 

1. For the Kuhn-Tucker sensitivity equations, most of the work in finding the 
parameter sensitivity is involved in the calculation (by finite differences) of the 
Hessian of the Lagrangian. However, after the first parameter sensitivity is 
determined the cost of evaluating successive sensitivity derivatives is reduced to 
(n+1) extra function evaluations. 

2. For the first order EDS algorithm, the work required to calculate the parameter 
sensitivity does not increase with problem size. However, this algorithm will 
not deliver 3u/dp and this algorithm may not be able to find the correct value for 
dx/8p. This will mean that df*/dp will also be inaccurate with this method. If 
the problem is fully constrained the accuracy of 9x/9p is better but the method 
may still provide inaccurate derivatives. 

3. For the second order EDS algorithm most of the work is in the calculation of the 
Hessian of the objective function and the Hessian of the constraints. The work 
involved for calculation of successive parameter sensitivities only requires 
approximately n+2 extra function evaluations. This algorithm requires the 
solution of a quadratic approximating problem for every new value of the 
parameter supplied by the user. 

4. For the RQP 1 algorithm (forward differencing) is the most efficient of the RQP 
methods proposed and seems to be much more efficient than the Kuhn-Tucker 
algorithm. The work required to calculate successive parameter derivatives is 
constant (n+1 function evaluations). This algorithm will perform well when B is 

a good approximation and the perturbation Ap is properly chosen. 

5. The number of function evaluations for the RQP 2 algorithm (forward 
differencing and 2 iterations of the RQOPT algorithm) grows linearly. The work 
required to calculate successive parameter derivatives is constant (2n+2 function 
evaluations). The work for calculating successive parameter sensitivities may be 
reduced because the Hessian approximation will improve after each parameter 
sensitivity derivative is approximated, which will eventually reduce the amount 
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of work required to solve the perturbed problem. 

6. For the RQP 3 algorithm (central differencing) the work involved grows linearly 
and the work for calculating successive parameter derivatives is constant (2n+2 
function evaluations). An indication of nonlinearity of the sensitivity derivatives 
can be indicated by checking for second derivatives of the functions as follows 


d 2f f + , 2f* + f- 
dp 2 Ap 2 


(3.46) 


This approximation of the second derivatives may not yield accurate results but it 
may be able to indicate that there is curvature present in the problem. Another 
advantage of using central differences occurs when the active set changes and 
directional derivatives can be approximated. 

7. The RQP algorithm with central differencing and 2 iterations of RQOPT is the 
most expensive of the proposed RQP algorithms. The work required to calculate 
successive parameter derivatives is constant (4n+2 function evaluations). The 
work for calculating successive parameter sensitivities will be reduced if we 
allow updating of the Hessian approximation during the RQP iterations, as less 
work will be required to solve the perturbed problem when the Hessian 
approximation is improved. 


The above discussion dealt with the number of required function evaluations to 
calculate the parameter sensitivities. We did not account for any of the other overhead such 
as solving the QP subproblem for the RQP method or solving a quadratic approximating 
problem for the second order EDS algorithm. 


The overhead associated with using the Kuhn-Tucker sensitivity equations is 
relatively small after the first parameter sensitivity is calculated, this is because if a 
factorization (i.e. LU) is used to solve the Kuhn-Tucker sensitivity equations then the 
amount of overhead becomes o(n) flops. The overhead for solving the RQP subproblems 
will also be realitively small if a good implementation of the RQP method is used (i.e. a 
proceedure propossed by Gill et. al. (1987) requires anly o(n) flops). The overhead for the 
first order EDS method will also be relatively small. However the overhead for the second 
order EDS method could be large depending on the problem. 

In summary, the RQP based methods are competitive with the existing methods. 

All variants of the RQP based method require approximately the same number of function 
evaluations for small problems (n<5), but considerably less for larger problems (n>5). 


3.4. POTENTIAL PROBT ,F.MS 


One of the main issues that needs to be investigated concerns the Hessian 
approximation: will the approximation converge in practice as predicted by the theory? If 
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convergence has not taken place then we need to investigate how to improve the Hessian 
approximation. Some modifications that can be made to obtain a more accurate Hessian 
approximation are discussed in Chapters 4 and 5. 

As with the estimation of any gradient by finite differences, the perturbation step 
size Ap and the nonlinearity of the problem will effect accuracy of the derivative 
approximation. Rules from Gill, Murray and Wright (1983) or Adelman, Haftka, and Iott 
(1986) can be investigated as a means to select the step size Ap. An automated selection 
proceedure for Ap should be investigated after the initial RQP sensitivity algorithm is 
tested. 


When using the forward difference option the choice of Ap is even more critical. If 
Ap is too small and the optimum of the problem is not known exactly then when the 
perturbed problem is solved we may only be seeing a better estimate of x* being found 
rather than an estimate of the solution of the perturbed problem. This will cause the 
derivative approximations to be inaccurate. If Ap is too large then we may only be 
obtaining trend information for the problem. 

All optimization programs incorporate some kind of convergence criteria that is 
based on the relative change in the design variables. This stopping criteria will effect the 
calculation of the sensitivity derivatives for all available methods, because there is a 
common assumption that the base point is a true optimum. The central difference 
approximation may be less sensitive to inexact solutions because the solution of the 
perturbed problems will be of a similar degree of accuracy. 

When solving the quadratic programming subproblem some type of convergence 
criteria is normally used. How small this tolerance is will effect how much work is needed 
to solve the subproblem (Nash 1985). During the early stages of the optimization it is not 
advisable to locate the exact solution of the QP subproblem as this may be too expensive. 
However once the program is in the region of a minimum the solution of the subproblem 
needs to be accurate. Therefore, we expect to use a tight convergence criteria for our QP 
solver during our reoptimizations. 


3.5. SUMMARY 

We have proposed a method and some variants based on the RQP method for 
estimating parameter sensitivities which provides sensitivity estimates nearly equivalent to 
the Kuhn-Tucker method. The method avoids the need for calculating second derivatives 
and its efficiency is competitive with current methods. The accuracy of the method 
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depends on two major pieces of information, the quality of the Hessian approximation 
provided by the RQP method, and the step size of the difference parameter used in the 
difference formula. Both these aspects of the method will be discussed in the following 
chapters. 
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4 . Implementation 


This chapter will discuss the implementation of the new parameter sensitivity 
method described in chapter 3. The program used as the basis for testing the new method 
was the RQOPT program which is an implementation of an active set RQP method 
(Beltracchi and Gabriele, 1987). The discussion begins with a discussion of the 
modifications made to RQOPT to perform the necessary calculations, and ends with a 
description of the software system developed to calculate parameter sensitivities. 

4.1. MODIFICATIONS TO RQOPT 

Most of the modifications to RQOPT were concentrated in one of the major areas of 
concern for the new sensitivity algorithm, the Hessian approximation. These modifications 
are discussed in subsections 4.1.1 and 4.1.2. The line search of RQOPT was also 
modified to yield a smoother convergence to the problem solution and this is discussed in 
subsection 4.1.3. The final modification discussed in subsection 4.1.4, provided the 
option of using a different variable metric update to yield a more accurate Hessian 
approximation 

4.1,1. Implementation of a Factorized BFS Variable Metric Update 

Variable metric updates have been successfully used for the past 20 years for 
unconstrained optimization and have been used successfully for approximately the past 10 
years for constrained optimization. Variable metric updates attempt to build an 
approximation to the Hessian matrix using only first order information, and solve for the 
search direction from the following equation 

s = B-!Vf (4.1) 

where B represents the approximation to the Hessian, Vf the gradient of the objective 
function, and s the search direction of the design variables. Variable metric updates have 
been provided in the literature for approximating either the inverse of the Hessian or the 
Hessian itself. 

Variable metric updates all have the same basic form. They begin with an 
approximation to the Hessian matrix, and then update the approximation by some rank one 
or rank two correction. The form of the update is normally 

Bnew = Bold + vv T + ww T (4.2) 

where v and w are calculated as some product of the old Hessian approximation, the last 
search direction, and the change in the gradient of the objective function. 


Several different forms of equation 4.2 have been proposed. The most popular 
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variable metric update has been the BFS (also known as the BFGS) which was proposed in 
1970 simultaneously by Broyden, Fletcher, Shanno and Goldfarb. The BFS update has 
been shown to be the best general purpose variable metric update. 

One of the problems associated with the BFS variable metric update is that it is 
effected by the problem scaling. Shanno and Phua (1978) have proposed a self scaling 
version of the BFS update. Its use in a RQP algorithm was investigated by Van der Hoek 
(1980). He found the self scaling variant with the second Oren-Spedicato (Oren 1974) 
switch seemed to perform the best with the particular RQP algorithm that he was using. 

In the mid 1970's several authors proposed updating the LDL t factors of the 
Hessian approximation with a procedure that could be used to stabilize the BFS update in 
terms of the numerical noise encountered in the calculation of the update. With the LDL T 
update we can be assured that the Hessian approximation remains positive definite, this will 
assure that the search directions that are generated from (4.1) are downhill. Additionally, 
finding the search direction from equation 4. 1 becomes a simple matrix calculation when 
using the LDL T update 

When variable metric updates are used for RQP methods it is normally preferred 
that the approximation of the Hessian of the Lagrangian be updated instead of its inverse. 
This is because solution of the QP subproblem requires the Hessian approximation. The 
BFS variable metric update is used by most of the successful implementations of the RQP 
method. 


The BFS update that was used in RQOPT is defined as 


B 


new 


_ Bo]d (zBqm )(zB 0 id) T | w w T 


z T B 0 idZ 


w T y 


where z and w are defined 




(4.3) 


z — x new ■ x old 


(4.4) 


y — ^ 7 xL( x new> v new> u new) - V x L(x 0 ld>Vnew> u new) 


(4.5) 


' 1 

© = - 0.8 z t Bz 
[z T Bz - z T y 


if z T y > 0.2 z T B z 
otherwise 


(4.6) 


w = © y + (1 - ©)Bz (4.7) 

Where the 0 term in equation 4.6 and 4.7 was defined by Powell (1977) to help maintain 
positive definiteness of the Hessian approximation, under normal operation 0 is equal to 
one. The Hessian approximation is guaranteed to be positive definite if z T w is greater than 
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zero. The Hessian approximation is not updated by RQOPT if z T w is less than zero. 

For this study, the LDL T update for the BFS variable metric (defined in equation 
4.3) as described by Gill and Murray (1978) was implemented (where z and w were 
calculated by equations 4.5 and 4.7). This update uses several matrix transformations to 
achieve a stable update. The actual update of the Hessian approximation is performed with 
a procedure described by Fletcher and Powell (1974) and extended by Gill, Murray, and 
Saunders (1975). 

In addition to the stability of this update relative to numerical noise, as discussed 
above, the LDL T update provides a convenient means for establishing a reset criteria for the 
Hessian approximation. The need for a reset of the Hessian approximation is discussed in 
the following section. 

4. 1 .2 Condition Number Reset 

Occasionally, due to numerical noise or a highly nonlinear problem, the Hessian 
approximation may become singular or indefinite. When this happens we can no longer be 
certain that the resulting search directions will satisfy the descent property that is assumed 
by the RQP. The only means to recover from this situation is to reset the approximation to 
some known positive definite matrix, which is generally the identity matrix. Early version 
of the BFS update were reset every n+1 iterations but this is a conservative approach that 
will sometimes erase good information and slow the convergence of the algorithm. The 
current thinking is to use a less conservative reset criteria that is based on a condition 
number estimate of the matrix with the hope that useful information built up in previous 
iterations is used for more iterations and should result in better convergence. 

The original version of RQOPT reset the Hessian approximation every time the 
active set changed or every n+1 iteration. A change in the active set results in a different 
QP subproblem to be solved and it was felt that the Hessian approximation would no 
longer be valid. Using this conservative reset criteria would prove unacceptable if we were 
using RQOPT to perform sensitivity analysis. With this reset criteria, we risk resetting the 
Hessian approximation just before the optimum is reached and would be left with only a 
few iterations of the method upon which to build an approximation. Thus we may have a 
very poor Hessian approximation when it comes time to perform the sensitivity analysis. 

The reset criteria adopted has been used successfully by several other algorithms 
(Powell 1985, Schittkowski 1983, Arora and Tseng 1987). The new reset criteria resets 
the Hessian approximation when the estimate of the condition number exceeds a fixed limit. 
This estimate can be found by computing 
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cond(H) est = 3 m22 - (4.8) 

a min 

where d m jn and d m ax are the smallest and largest values of the D matrix in the LDL T 
factorization. 

Using this reset criteria has led to a more stable update yielding faster convergence 
for the RQOPT program and more accurate estimates of the Hessian of the Lagrangian. 

4.1.3 Calculation of the Lagrange Multiplier Estimates 

The Lagrange multiplier estimates are an integral part of building the Hessian 
approximation. The value of the Lagrange multiplier estimates are used as inputs to the 
variable metric update to approximate the Hessian of the Lagrangian function. 

The original version of RQOPT calculated the Lagrange multiplier estimates as the 
Lagrange multipliers of the constraints in the QP subproblem. This value of the Lagrange 
multiplier estimate is a valid estimate of the true multipliers when a step of a = 1 is used in 
the line search (Gill and Murray 1979). When this occurs, the estimates should converge 
to the true Lagrange multipliers as the problem converges. 

A problem can arise, however, in the first few iterations of RQOPT. At the 
beginning of a search it is possible that a Lagrange multiplier estimates produced by the QP 
subproblem will be several orders of magnitude larger than true value of the Lagrange 
multiplier. If the line search then makes a small step (a « 1), the large value of the 
Lagrange multiplier estimate may bias the updating of the Hessian approximation in such a 
way that new approximation only sees the constraint associated with the large Lagrange 
multiplier. It may then take several iterations before the Hessian approximation is 
corrected. 

RQOPT was modified to use the following linear interpolation to update the value 
of the Lagrange multiplier estimates after the line search is completed 

u new = uold + Oc(llqp - Uold) (4.9) 

When a step length of a =1 is used in the line search (equation 3.4) then formula 4.9 
updates the Lagrange multiplier estimates to be the estimates delivered by the QP 
subproblem. This update was also used by Schittkowski (1983). 

The procedure for updating the Lagrange multiplier estimates helped yield a 
smoother convergence of the Hessian approximation, because we were able to more 
accurately represent the Lagrangian function when we were performing the approximation 
updates. 
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4.1.4. SRI update 


The SRI update is a variable metric update that does not require exact line searches 
for quadratic convergence, where as the BFS update requires exact line searches for 
quadratic convergence. Because the RQP method seldom performs exact line searches, it 
was felt the SRI update may be able to obtain a better approximation of the Hessian of the 
Lagrangian. 

A table describing the differences between the BFS and SRI update is presented 

below 

Update Advantages Disadvantages 

BFS Self Correcting Requires exact line searches 

Stable (maintains positive definiteness) 

Has a good performance history 

SR 1 Does not require exact line searches update may be undefined and it is not 

guaranteed to maintain positive definiteness 
of the Hessian Approximation. There is not 
a lot of literature on the performance of this 
update. 

Table 4.1 A comparison of the BFS and SRI variable metric updates 

The stability of the BFS variable metric update has led to its use in almost all RQP 
implementations. However Cha and Mayne (1987) report that they have tested the SRI 
update and found exact convergence of Hessian approximations for quadratic functions. 
Although the SRI update lacks the stability of the BFS update, we were interested in 
comparing the performance of the 2 updates in terms of the Hessian convergence. If the 
SRI update delivers better Hessian approximations than the BFS update then we will have 
to further investigate methods to stabilize the SRI update. 

The SRI update is defined as follows 

b -=^ | (B, ^)' 2>T (410) 

where y and z are obtained from equation 4.4 and 4.5. This update is undefined when the 

denominator is equal to zero. The SRI update may be undefined even for positive definite 
quadratic problems. This problem was addressed by Brayton and Cullem (1979), Cullem 
and Brayton (1979). 

The symmetric rank one (SRI) update was implemented in both a factored (LDL t ) 
and unfactored form. In our implementation if the absolute value of the denominator (in 
equation 4.10) is less than some small number we use the BFS update which is described 
in section 4.1.1. 
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Even though the SRI update may be undefined, it has the very nice property of not 
requiring exact line searches. This is important because in the RQP method we do not 
perform exact line searches, and the BFS variable metric method assumes exact line 
searches. Powell (1986) clearly demonstrates the detrimental effect of inexact line searches 
on the BFS method. The performance of the SRI update for solving quadratic problems is 
such that after n updates (providing that all updates are defined) the Hessian approximation 
will have converged to the true Hessian. Thus we may obtain a better convergence of the 
approximation of the Hessian of the Lagrangian if we are able to use the SRI update. 

Some preliminary results were obtained comparing the BFS and SRI variable 
metric updates and these are discussed in section 5.3. 

4.2 THE CREATION OF A SYSTEM TO AUTOMATICALLY CALCULATE 
PARAMETER SENSITIVITIES 


In this section we provide a brief overview of the software system created for 
studying parameter sensitivities. The software system is made up of three major pieces: a 
problem preparation package RQCRE, the RQP algorithm using the modifications 
described in the preceding section, RQOPT, and an interactive program RQSEN, that acts 
as a post processor/sensitivity analysis module for the RQOPT program. The RQOPT 
program was an existing program and has been documented previously (Beltracchi and 
Gabriele, 1986). The RQCRE and RQSEN programs were created for this study and will 
be briefly described in the following paragraphs. A more detailed discussion of these 
systems is provided in the appendix. 

4.2. 1 The RQCRE Support System 

The RQCRE program is set up to be used as an interactive tool for use with the 
RQSEN system. The purpose of the RQCRE program is to remove the chance of errors in 
the problem formulation. The RQSEN program requires approximately 30 arrays to be 
dimensioned which are automatically dimensions by RQCRE. The RQCRE program also 
automatically writes the calling program and data files required by the RQSEN system. 

The RQCRE program requires the user to provide basic information about the 
problem such as the number of variables, number of equality constraints, number of 
inequality constraint, and number of parameters that will be studied. 

The RQCRE system then produces a main calling program, a shell of the function 
subprogram 1 used to define the objective function and the constraints, and a data file used 
for input into the RQSEN system (sample output is provided in the appendix). The 

1 The RQCRE program is not designed to allow the user to enter definitions of the objective function or 
constraints, these definitions must be entered manually into the code that was generated by RQCRE. 
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RQCRE program also sets up the default values for the algorithm parameters used by 
RQOPT. 

4.2.2 The ROSEN program 

The RQSEN system was set up as a pre and post processor for the RQOPT 
program. The RQSEN system was set up to be an interactive user friendly program for 
performing the following basic functions; 

1. The system can be used to solve optimal design problems 

2. The system can be used to calculate parameter sensitivities 

3. The system can be used to conduct studies of large variations in problem 
parameters 

4. The system is also set up to create sensitivity plots of that can be used to 
perform trade off studies. 

A sample session with the RQSEN system illustrating these options is presented in the 
appendix. 

The RQSEN system requires a calling program and a function subprogram 
(defining the objective function and the constraints) to be written in FORTRAN 2 . The 
RQSEN system also requires the user to define a data file that contains the algorithm 
parameters, and the initial values of the design variables and design parameters. The user 
can then direct the RQSEN program to study the sensitivities of only certain parameters. 

The RQSEN program first produces optimum designs. Once the problem has been 
optimized the RQSEN system can be used to produce parameter sensitivity derivatives, 
which can then be used to study the effect on the optimum of large variations in the 
parameters. The RQSEN system is also set up so that an external graphing program can be 
used to create plots of the optimum sensitivities for large variations in the parameters can be 
studied. A typical plot is presented in figure 

2 The RQCRE system can be used as an aid in creating the calling program and function subprogram. 
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p(3) with p(l)=3, p(2)=l 

Figure 4. 1 A plot of the Sensitivity of the Optimum of test problem 1 to p(3) 


Plots similar to this one can also be generated for the design variables, Lagrange 
multipliers and values of the constraints. These plots can then be used to assess the 
characteristics of the problem (such as nonlinearity and changes in the active set). Using 
these plots to assess the characteristics of the problem will be discussed in the results part 
of chapter 5. Plots similar to figure 4.1 are presented in the appendix for problems in the 
test set. 
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5 . 


Numerical Experiments 


This chapter describes the numerical experiments that have been conducted to date 
on the new sensitivity method. We begin by discussing the initial test set used and any 
special features of the selected problems. Next, we discuss testing that has been performed 
comparing the accuracy of the known Hessian to the approximations obtained, which 
includes comparisons of the BFS and SRI updates. In the third section, the accuracy of 
the sensitivity derivatives obtained with the new sensitivity algorithm is assessed against 
the known results. This section also compares the effect of choosing a central or forward 
difference formula and the effect of the step size Ap. The final section presents some 
conclusions drawn from this initial testing. 

5.1. INITIAL TEST SET 

A two phase testing program has been formulated for studying the effectiveness of 
the new method for estimating parameter sensitivity. The first phase was to develop a set 
of test problems for which the parameter sensitivities could be exactly determined using the 
Kuhn-Tucker equations. This required that any second order information needed could be 
determined analytically. Choosing problems of this type would allow a direct comparison 
of the sensitivity results produced by the new method with the exact sensitivities and also 
allow the comparison of the BFS and SRI Hessian approximations. From this study we 
hope to develop some insight into several questions concerning the algorithm such as: 
proper choices for algorithm parameters (i.e. the proper size of Ap), what is the most 
reliable Hessian approximation, how close does the Hessian approximation have to be to 
achieve good results, does updating the Hessian approximation during the sensitivity 
analysis significantly improve the estimate, and which of the variants (forward/central 
difference approximations with one or two iterations of RQOPT) described in chapter 3 
provides the most consistent results. 

The second phase of the testing would consist of testing the algorithm against a set 
of engineering problems where second order information would not be available. Here the 
results obtained from the sensitivity algorithm would be compared to actual reoptimization 
results to assess its accuracy. In the time allotted for this study, only the first phase of 
testing has been completed and is reported on here. 

The problems making up the initial test set are presented in the appendix of this 
report. We have experimented so far with 4 test problems that have a total of 12 
parameters. The problems possess both linear and nonlinear behavior. We expect to 
expand this test set in the near future. Plots of the optimum sensitivity for selected 
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problems and parameters are also presented in the appendix. 


5.2. CONVERGENCE OF THE HESSIAN APPROXIMATION 

The derivation given in section 3.2.1 showed that equivalence of the new method 
with the Kuhn-Tucker method depends on the accuracy of the Hessian approximation 
obtained from the RQP method. Using this initial test set we hope to observe how closely 
the Hessian approximation comes to the exact Hessian and draw some initial conclusions 
on its importance to the accuracy of the results. 

A measure of the closeness of the Hessian approximation to the true Hessian can be 
defined using the Frobenius norm as 

Eh = II H - Happrox Hf (5.1) 

This measure has been used in the past to compare the convergence of different variable 
metric updates (Dennis and Schnable 1983). 

For test problem 1 the true Hessian of the Lagrangian is 

H-P 64 0 1 

H_ L 0 2.6 J 

From the RQOPT program we obtained the following Hessian approximation with the BFS 
update 

„ _r 1.50017 -.5403101 
H BFS~L_ 540310 2.34388 J 

which gave us a ehbfs = 1-396. 

Using the SRI update from the same starting point, we obtained the following 
Hessian approximation 

„ r 2.63194 -.002930 1 
Hsri_ L-.002930 2.61374 J 

with gives a Ehsri = 0.0164. This represents a large improvement in the closeness of the 
Hessian approximation. However, even though the Hessian approximation for the SRI 
update is much better than the Hessian approximation for the BFS update the problems 
were solved in the same number of iterations (and function/constraint evaluations) of 
RQOPT. 

The results given above were obtained with a value of 8 =1.1. The 8 parameter 
controls the size of the active set during the course of an optimization; a large 8 will cause 
more near active constraints to be considered as part of the active set, a small value of 8 will 
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allow only truly active constraints to be considered. Having the proper set of active 
constraints identified early in the optimization could effect the accuracy of the Hessian 
approximation. To test this, a larger value of 8 (8 =10.1) was chosen and the problem 
resolved obtaining the following Hessian approximations 

„ T2.329 .3227 1 

Hbfs_ L.3227 2.264 J 

„ JT2.6394 .00093 l 

h SR1-|_.00093 2.5990 J 

the values of ehbfs = 0.6019 and Ehsri = 0.00174 were obtained. These improved 
Hessian approximations result because RQOPT was able to identify the correct active set of 
constraints sooner. With the large value of 8 RQOPT required the same number of 
iterations to solve the problem, but required more constraint evaluations. 


Another implementation issue that needs investigation concerns whether the 
Hessian approximation obtained from the optimization should be further updated during the 
reoptimizations performed to estimate the sensitivities. To study the effect of allowing 
Hessian updates during the reoptimization, the sensitivity with respect to parameter 3 in 
problem 1 was estimated with this option enabled. The Hessian approximation that was 
used at the start of the sensitivity analysis is the Hessian approximation that was obtained 
with the BFS update and 8 = 1.1. After estimating the sensitivity, we obtained the 
following Hessian approximation 

„ _r2.63975 .00013 1 

h BFS-[. 00013 2.59957 J 

with Ehbfs = 0.00053. This indicates that there is a possibility for improving the Hessian 
approximation if we allow updating during the sensitivity analysis. 


Tests for problem 2 were also performed, whose true Hessian of the Lagrangian is 
given by 


H= 


r5 i 1 I 
1 5 1 
Ll 1 5 J 


Using the starting point provided in the problem description, we obtain the following value 
of the Hessian approximation (from RQOPT) when we use the BFS update 


Hbfs= 


4.4976 1.9976 0.49763-1 
1.9976 2.9976 1.9976 
.0.49763 1.9976 4.4976 . 


with ehbfs = 3.000. When we use the SRI update we obtain the following value of the 


Hessian approximation 
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Hsri= 


[4.5 2.0 0.5 ~\ 
2.0 3.0 2.0 
l0.5 2.0 4.5. 


£Hsri = 3.0. If we allow the Hessian matrix to be updated while estimating the sensitivity 
of pi with a Ap = 0.0001 we obtain 


Hbfs= 


[4.9448 1.0070 1.1749 1 
1.0070 4.9964 0.9665 
Ll. 1749 0.9665 4.3990 J 


Hsr1 = 


r511i 
1 5 1 

Ll 1 5 J 


with ehbfs = 0.6541 and with ehsri = 0.00. This represents a significant improvement 
of the Hessian approximations, particularly when the SRI update is used. 


If we calculate the sensitivity of p 2 and use the SRI update we also obtain exact 
convergence of the Hessian approximation. However if we use the BFS update we do not 
obtain exact convergence but an improvement similar to that of the first problem is 
achieved. 


For Test problem 3 the Hessian of the Lagrangian is 


H= 


[12 0 0 0 "I 
0 8 0 0 
0 0 10 0 
Lo 0 0 4 J 


If we use the starting point that was provided in the problem description, the approximation 
to the Hessian of the Lagrangian (form RQOPT) using the BFS update is 


r 9.785 -0.4657 -2.502 -0.9879 1 


Hbfs= 


-0.4657 7.7556 -0.5500 0.0174 
-2.502 -0.5500 4.062 0.7464 

L-0.9879 0.0174 0.7464 2.1579 J 


with Ehbfs = 7.76. If we use the SRI update to solve the problem then we obtain the 
following approximation to the Hessian 


Hsri= 


[11.9744 -0.02493 -0.04194 0.02212 i 
-0.02493 7.9792 -0.03037 0.01095 
-0.04194-0.03037 9.9679 0.00222 
L 0.02212 0.01095 0.00222 4.01399 J 


with a Ehsri = 0.130. This represents a major improvement in the Frobenius norm. 


For Test problem 4 the Hessian of the Lagrangian is 


r6.72 

-4.0 

-2.0 

6.4 

-2.0 

-4.0 < 

O 

o 

ON 

-1.2 

-6.2 

6.4 

-2.0 

-1.2 

4.4 

-1.2 

-2.0 

6.4 

-6.2 

-1.2 < 

L3418 

-4.0 

L-2.0 

6.4 

-2.0 

-4.0 

6.2688 
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using the starting point that was defined in the appendix, RQOPT with the BFS update 
yields the following Hessian approximation 


r 6.280 -3.963 
-3.963 8.052 


HbfsH 


-1.417 -0.561 
6.458 -6.247 
L- 1.924 5.775 


-1.417 6.458 -1.924-1 
-0.561 -6.247 5.775 
1.548 -0.932 -1.449 
-0.932 9.3465 -4.024 
-1.449 -4.024 5.974 J 


with £hbfs = 3.6226. 


When we attempted to use the SRI update, the Hessian approximation became 
nearly singular after 5 iterations and the Hessian approximation was automatically reset to 
the identity matrix by RQOPT. RQOPT delivered the following Hessian approximation 


r 4.661 1 -4.8848 0.1290 5.2854 -3.5329 
-4.8848 7.5178 -.1721 -7.0517 4.7140 
H S ri= 0.1290 -0.1721 1.0046 .1862 -.1245 
5.2854 -7.0517 0.1862 8.6328 -5.0996 
L- 3.5329 4.7140 -.1245 -5.0996 4.4094 


with £hsri - 9.307. The inaccuracy of this Hessian approximation results because a total 
of only 7 iterations were needed to solve the problem, and a reset occurred after the fifth 
iteration. Therefore, only 2 iterations could be used to build the Hessian approximation. 
In the near future we will investigate why the Hessian approximation became nearly 
singular after the 5 lh iteration. 


A summary of the results of this section are presented in the Table 5. 1 where £o 
represents the error between the true Hessian and the identity matrix used at the outset of 
the optimization. Using the BFS update we see that we were not able to converge to the 
exact Hessian but the inaccuracies do not seem to be large. As mentioned before, this may 
be due to RQOPT not using exact line searches which the BFS method assumes. Allowing 
updating of the Hessian approximation during the sensitivity analysis seems to improve the 
estimate of the Hessian of the Lagrangian. 


Using the SRI update we were able to obtain better estimates of the Hessian of the 
Lagrangian for both problem 1 and 3. For problem 2 the Hessian of the Lagrangian that 
was produced by the SRI update had converged in a projected or reduced sense. The 
inaccuracies in problem 4 are due to a near singular point which is discussed above. 


Problem 

£0 

£BFS 

£SRi 

e BFS with updating 

1 

2.291 

1.396 

0.061 

0.0005 

2 

7.348 

3.0000 

3.0 

0.654 

3 

15.93 

7.76 

0.130 

5.186 

4 

23.276 

3.6226 

9.307 

1.698 


Table 5.1 A comparison of the Frobenius norms of the Hessian approximations 
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5.3. RESULTS 


This section presents a comparison of the sensitivity derivatives calculated by the 
new method with the known sensitivities for the problems in the initial test set. We also 
present a means that can be used to compare the accuracy of the sensitivity derivatives 
graphically. 

The measure for accuracy that will be used was also used by Sandgren (1977). 
Sandgren compared the closeness of the optimum design point generated to known 
optimum point, and the closeness of the value of the known optimum objective function 
value to the generated value of the optimum plus a penalty for any violated constraints. 
Sandgren defined the following measures 

e f for f(x * 

L\BS[f(x)] for f(x* 


) * 0 
) = 0 


(5.2) 


where f(x*) is the true value of the optimum and f(x) is the value returned by the 
algorithm. The total error is calculated as 

nineq neq 

e t = £f+ I<gj>+ I(hi) (5.3) 

j=l i=l 

where < a > = (0, if a > 0 1 -a if a < 0). The e t measure is used because it does not bias any 
constraints. 


The relative error in the x vector is defined as 



in equation 5.4 if xj* is equal to zero then the relative error in xi is defined as the value of 

Xi. 


We will define the relative error in the gradient (df*/dp) of the objective function as 

follows 
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r 


ABS 


e df*/dp 




■df* df* - 
■3p st ' dp 
dP 
dp 


df * 

f « r -hr* 0 

dp 


(5.5) 


Hf~] ford ^=° 

We will define the relative error in dx*/3p and 3u*/9p in the same manner as e x and denote 
these values as £dx*/dp and £0 u */dp respectively. Eight digits of accuracy were maintained 
in calculating the relative errors. 


The optimal sensitivities for the test problems were calculated using the Kuhn- 
Tucker method with exact derivatives. Once the optimal sensitivities for the problems were 
known, experiments were conducted using RQSEN on the initial test set. Both the forward 
difference and central difference variants of the RQP sensitivity algorithm were tested with 
large and small values of perturbation for the parameters. For all cases, RQOPT 1 was 
allowed to perform two iterations to optimize the perturbed problem. However, there were 
some instances where RQOPT required only one iteration to meet the convergence criteria. 
A spreadsheet was used to automate the calculation of the relative errors in the derivatives 
using the formulas given above. Summary tables showing the relative errors in the 
calculation of the derivatives will be presented for each of the problems. 


Plots of the optimal sensitivities for large variations in the parameters were also 
created for all of the parameter sensitivities that were studied. The interesting plots will be 
included in the appendix of this report. These plots can be used to help asses the 
nonlinearity in the sensitivity derivatives, and to help to understand the effect of changes in 
the active set. 


The rest of this section presents tables and figures showing the relative accuracy of 
the sensitivity derivatives. A brief discussion of the results for each problem is offered. 

5.3.1 Problem 1 

Problem 1 possesses three parameters for study. Sensitivities of the objective 
function, design variable, and Lagrange multipliers for each parameter were estimated 
using the four variations of the basic algorithm. The results are compared against the exact 
sensitivities in Tables 5.2 - 5.7. In most cases, the estimated sensitivities agree with the 
known sensitivities with few exceptions. As might be expected, the central difference 
approximations in all cases provides better estimates than the forward difference 
approximations. No strong conclusions with respect to the choice of Ap can be drawn 
from this problem. For parameter 1, both sizes of Ap provide exact sensitivities. For 

*The gradients of the objective function and constraints were calculated using central and forward difference 
approximations. 


51 




parameter 2, the larger value of Ap provides better results while for parameter 3, the 
smaller value of Ap provides the best results. A review of the sensitivity plot for this 
parameter (figure A.2) shows that the sensitivities for this parameter are nonlinear. 

In conclusion, for this problem, using a central difference approximation with either 
step size for Ap resulted in no significant errors in the sensitivity estimates. 


Kuhn Tucker Central Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

1.0000000 

1.00000000 

0.00E+00 

1.00000000 

0.00E+00 

dxi/dp 

dx^/dp 

0.0000000 

0.0000000 

0.000000 

0.000000 

0.000E+00 

0.000E+00 

0.000000 

0.000000 

0.000E+00 

0.000E+00 

Ex 



0.00E+00 


0.00E+00 

dui/dp 

duVdp 

-0.2000000 

0.4000000 

-0.200000 

0.400000 

0.000E+00 

0.000E+00 

-0.200000 

0.400000 

0.000E+00 

0.000E+00 

Eu 



0.00E+00 


0.00E+00 


Table 5.2 Central Difference Approximations to the Parameter Sensitivities for problem 1 
parameter 1 


Kuhn Tucker Forward Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

1.0000000 

1.06000000 

-6.00E-02 

1.00300000 

-3.00E-03 

dxi/dp 

dx^/dp 

0.0000000 

0.0000000 

0.000000 

0.000000 

0.000E+00 

0.000E+00 

0.000000 

0.000000 

0.000E+00 

0.000E+00 

Ex 



0.00E+00 


0.00E+00 

dui/dp 

du^dp 

-0.2000000 

0.4000000 

-0.2000301 

0.39984547 

-1.503E-04 

3.863E-04 

-0.2006318 

0.39693215 

-3.159E-03 

7.670E-03 

Eu 



4.15E-04 


8.29E-03 


Table 5.3 Forward Difference Approximations to the Parameter Sensitivities for problem 1 
parameter 1 
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Kuhn Tucker Central Difference Approximations 


Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp -0.3000000 

-0.30000000 

0.00E+00 

-0.30000000 

0.00E+00 

dxj/dp 0.1000000 
d X2 /dp -0.2000000 

0.10000011 

-0.20000030 

-1.100E-06 

3.000E-06 

0.10000359 

-0.20000320 

-3.590E-05 

-1.600E-05 

£x 


3.20E-06 


3.93E-05 

dui/dp -0.1304000 
d U2 /dp 0.0008000 

-0.13040090 

0.00079973 

-6.902E-06 

3.363E-04 

-0.13040262 

0.00080051 

-2.009E-05 

-6.325E-04 

£u 


3.36E-04 


6.33E-04 


Table 5.4 Central Difference Approximations to the Parameter Sensitivities for problem 1 
parameter 2 


Kuhn Tucker Forward Difference Approximations 


Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp -0.3000000 

-0.30000000 

0.00E+00 

-0.30000000 

0.00E+00 

dxj/dp 0.1000000 
d X2 /dp -0.2000000 

0.10004811 

-0.20017630 

-4.811E-04 

-8.815E-04 

0.10000681 

-0.20000665 

-6.810E-05 

-3.325E-05 

e x 


1.00E-03 


7.58E-05 

d Ul /dp -0.1304000 
d U2 /dp 0.0008000 

-0.13061831 

0.00114555 

-1.674E-03 

-4.319E-01 

-0.12851547 

0.00999551 

1.445E-02 

-1.149E+01 

£u 


4.32E-01 


1.15E+01 


Table 5.5 Forward Difference Approximations to the Parameter Sensitivities for problem 1 
parameter 2 
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Kuhn Tucker Central Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

-0.4000000 

-0.40000000 

0.00E+00 

-0.40000000 

0.00E+00 

dxi/dp 

dx2/dp 

1.2000000 

0.6000000 

1.19995460 

0.60007142 

3.783E-05 

-5.952E-05 

1.19999990 

0.60000018 

1.000E-06 

-3.000E-07 

Ex 



7.05E-05 


1.04E-06 

dui/dp 

du2/dp 

0.4048000 

-0.5896000 

0.40501856 

-0.58972201 

-5.399E-04 

-2.069E-04 

0.40480055 

-0.58960030 

-1.359E-06 

-5.088E-07 

£u 



5.78E-04 


1.45E-06 


Table 5.6 Central Difference Approximations to the Parameter Sensitivities for problem 1 
parameter 3 


Kuhn Tucker Forward Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

-0.4000000 

-0.36400000 

9.00E-02 

0.39820000 

9.00E-02 

dxj/dp 

dx^/dp 

1.2000000 

0.6000000 

1.20029290 

0.59484082 

-2.441E-04 

8.599E-03 

1.20001670 

0.59973857 

-1.392E-05 

4.357E-04 

Ex 



8.60E-03 


4.36E-04 

dui/dp 

du^/dp 

0.4048000 

-0.5896000 

0.39514199 

-0.57797893 

2.386E-02 

1.971E-02 

0.40367242 

-0.59207284 

2.786E-03 

-4.194E-03 

£u 



3.09E-02 


5.03E-03 


Table 5.7 Forward Difference Approximations to the Parameter Sensitivities for problem 1 
parameter 3 

5.3.2 Problem 2 


Tables 5.8 - 5. 1 1 present results obtained for the two parameters of problem 2. 
There is some significant disagreement between the known and estimated sensitivities for 
parameter 1 in each of the variations tested. The inaccuracies seem to occur due to the 
Hessian approximation not converging. The cause of this is likely due to RQOPT not 
using exact line searches, which the BFS variable metric update assumes. In this case, 
possibly the SRI update would produce better results. The results for parameter 2 are 
excellent for all variations. 

The major conclusion to draw from this problem is that the convergence of the 
Hessian approximation can be a critical factor in the success of the new method. 
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Kuhn Tucker 
Method 

Central Difference Approximations 
AP =2% relative error AP=0. 1 % 

relative error 

df/dp -7.0000000 

-7.00000000 

-1.43E-11 

-7.00000000 

-1.43E-11 

dxi/dp 2.0888889 
d X2 /dp -2.1666667 
dx 3 /dp -0.9166667 

2.1593080 

-2.1459360 

-1.0143604 

-3.371E-02 

9.568E-03 

-1.066E-01 

2.1598250 

-2.1459879 

-1.0139152 

-3.396E-02 

9.544E-03 

-1.061E-01 

£x 


1.12E-01 


1.12E-01 

dui/dp -4.6666667 

-4.5219200 

3.102E-02 

-4.5311680 

2.904E-02 

£u 


3.10E-02 


2.90E-02 

dg2/dp -6.00000000 

-6.1756450 

-2.927E-02 

-6.1738965 

-2.898E-02 

e g 2.93E-02 2.90E-028 

Table 5.8 Central Difference Approximations to the Parameter Sensitivities for problem 2 
parameter 1 

Kuhn Tucker 
Method 

Forward Difference Approximations 
AP=2. % relative error AP=0. 1 % 

relative error 

df/dp -7.0000000 

-7.00000000 

-1.43E-11 

-7.00000000 

- 1.43 E- 11 

dxi/dp 2.0888889 
dx^dp -2.1666667 
dx 3 /dp -0.9166667 

2.2087420 

-2.1161414 

-1.1367750 

-5.738E-02 

2.332E-02 

-2.401E-01 

2.2379818 

-2.1100841 

-1.1301356 

-5.738E-02 

2.332E-02 

-2.401E-01 

Ex 


2.48E-01 


2.48E-01 

dui/dp -4.6666667 

-4.4230730 

5.220E-02 

-4.3803787 

5.220E-02 

Eu 


5.22E-02 


5.22E-02 

dg2/dp -6.00000000 

-6.43386710 

-7.231E-02 

-6.3725934 

-7.231E-02 


e g 


7.23E-02 


7.23E-02 


Table 5.9 Forward Difference Approximations to the Parameter Sensitivities for problem 2 
parameter 1 
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Kuhn Tucker 

Central 

Difference Approximations 


Method 

AP =2% 

relative error 

AP=0.1% 

relative error 

df/dp 

12.00000000 

12.00000000 

0.00E+00 

12.00000000 

0.00E+00 

dxi/dp 

dx^dp 

dx3/dp 

0.33333333 

0.33333333 

0.33333333 

0.3333333 

0.3333333 

0.3333333 

1.000E-08 

1.000E-08 

1.000E-08 

0.3333333 

0.3333333 

0.3333333 

1.000E-08 

1.000E-08 

1.000E-08 

Ex 



1.73E-08 


1.73E-08 

dui/dp 

2.33333333 

2.3333333 

1.429E-08 

2.3333335 

-7.143E-08 

e u 



1.4 3 E -08 


7.14E-08 

dg2/dp 

2.00000000 

2.0000000 

O.OOOE+OO 

2.0000000 

0.000E+00 

e g 



O.OOE+OO 


0.00E+00 


Table 5.10 Central Difference Approximations to the Parameter Sensitivities for problem 2 
parameter 2 


Kuhn Tucker Forward Difference Approximations 



Method 

AP=2.% 

relative error 

AP=0.1% 

relative error 

df/dp 

12.00000000 

12.00000000 

0.00E+00 

12.00000000 

0.00E+00 

dxj/dp 

dx^/dp 

dx3/dp 

0.33333333 

0.33333333 

0.33333333 

0.3333333 

0.3333333 

0.3333333 

1.000E-08 

1.000E-08 

1.000E-08 

0.3333333 

0.3333336 

0.3333333 

1.000E-08 

1.000E-08 

1.000E-08 

£x 



1.7 3 E -08 


1.73E-08 

dui/dp 

2.33333333 

2.3333333 

1.429E-08 

2.3333335 

1.429E-08 

e u 



1.43E-08 


1.43E-08 

dg2/dp 

2.00000000 

2.0000000 

0.000E+00 

2.0000000 

0.000E+00 

£ g 



0.00E+00 


0.00E+00 


Table 5.1 1 Forward Difference Approximations to the Parameter Sensitivities for problem 2 
parameter 2 

5.3.3 Problem 3 


Problem 3 possesses 3 parameters for study, and the results are presented in Tables 
5. 12 - 5. 17. Again, the central difference approximations produce the better results, with 
the small step perturbation better than the large step perturbation for the first two 
parameters. The accuracy of the estimates are good in comparison with the known 
sensitivities. 
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Kuhn Tucker 

Central 

Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

-8.0000000 

-8.0000001 

-8.75E-09 

-8.0000001 

-8.75E-09 

dxi/dp 

dx2/dp 

dx3/dp 

dx^dp 

-1.1471984 

-0.3342830 

-0.2279022 

-3.5403608 

-1.1551824 

-0.3360517 

-0.2189132 

-3.5312020 

-6.960E-03 

-5.291E-03 

3.944E-02 

2.587E-03 

-1.1485459 

-0.3355475 

-0.2265911 

-3.5390270 

-1.175E-03 

-3.783E-03 

5.753E-03 

3.767E-04 

£x 



4.05E-02 


6.99E-03 

dui/dp -67.3428300 
du 3 /dp 55.4605800 

-67.5423370 

55.6363930 

-2.963E-03 

-3.170E-03 

-67.3399900 

55.5155000 

4.217E-05 

-9.903E-04 

e u 



4.34E-03 


9.91E-04 

dg2/dp 

-1.6600260 

-1.6579776 

1.234E-03 

-1.6595000 

3.169E-04 

e g 



1.23E-03 


3.17E-04 


Table 5.12 Central Difference Approximations to the Parameter Sensitivities for problem 3 
parameter 1 



Kuhn Tucker 

Forward Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

-8.0000000 

-8.0000001 

-8.75E-09 

-8.0000001 

-8.75E-09 

dxj/dp 

dx^dp 

dx 3 /dp 

dxVdp 

-1.1471984 

-0.3342830 

-0.2279022 

-3.5403608 

-1.0434709 

-0.3298696 

-0.2847166 

-3.5084682 

9.042E-02 

1.320E-02 

-2.493E-01 

9.008E-03 

-1.1427923 

-0.3378248 

-0.2292447 

-3.5376210 

3.841E-03 

-1.060E-02 

-5.891E-03 

7.739E-04 

EX 



2.66E-01 


1.27E-02 

dui/dp -67.3428300 
du 3 /dp 55.4605800 

-62.2695100 

52.1090250 

7.534E-02 

6.043E-02 

-66.8998810 

55.1460460 

6.578E-03 

5.671E-03 

Eu 



9.66E-02 


8.68E-03 

dg2/dp 

-1.6600260 

-1.6647592 

-2.851E-03 

-1.6589239 

6.639E-04 

e g 



2.85E-03 


6.64E-04 


Table 5.13 Forward Difference Approximations to the Parameter Sensitivities for problem 3 
parameter 1 
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Kuhn Tucker Central Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

0.0000000 

0.0000000 

0.00E+00 

0.0000000 

0.00E+00 

dxi/dp 

dx2/dp 

dx3/dp 

dxVdp 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.000E+00 

0.000E+00 

0.000E+00 

0.000E+00 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.000E+00 

0.000E+00 

0.000E+00 

0.000E+00 

ex 



0.00E+00 


0.00E+00 

dui/dp 

du3/dp 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.000E+00 

0.000E+00 

0.0000000 

0.0000000 

0.000E+00 

0.000E+00 

£u 



0.00E+00 


0.00E+00 

dg2/dp 

1.00000000 

1.0000000 

0.000E+00 

1.0000000 

0.000E+00 

e g 



0.00E+00 


O.OOE+OO 


Table 5.14 Central Difference Approximations to the Parameter Sensitivities for problem 3 
parameter 2 


Kuhn Tucker Forward Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

0.0000000 

-2.200E-13 

2.200E-13 

-1.660E-12 

1.660E-12 

dxi/dp 

dx^dp 

dx3/dp 

dx^dp 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

4.150E-05 

-3.300E-05 

-2.020E-05 

-3.090E-05 

-4.150E-05 

3.300E-05 

2.020E-05 

3.090E-05 

8.316E-04 

-6.601E-04 

-4.050E-04 

-6.180E-04 

-8.316E-04 

6.601E-04 

4.050E-04 

6.180E-04 

£x 



6.46E-05 


1.29E-03 

dui/dp 

du3/dp 

0.0000000 

0.0000000 

8.840E-03 

-1.270E-02 

-8.840E-03 

1.270E-02 

1.769E-01 

-2.553E-01 

-1.769E-01 

2.553E-01 

£u 



1.55E-02 


3.11E-01 

dg2/dp 

1.00000000 

1.0000100 

-1.000E-05 

1.0020100 

2.010E-03 

e g 



1.00E-05 


2.01E-03 


Table 5.15 Forward Difference Approximations to the Parameter Sensitivities for problem 3 
parameter 2 
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Kuhn Tucker Central Difference Approximations 


Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp -10.0000000 

-9.99999986 

1.4E-08 

-9.99999986 

1.4E-08 

dxi/dp 1.3057920 
dx2/dp 0.5460588 
dx 3 /dp -1.0541310 
dxVdp -2.3741690 

1.3055639 

0.5465549 

-1.0523005 

-2.3704118 

1.747E-04 

-9.085E-04 

1.737E-03 

1.583E-03 

1.3079005 

0.5481210 

-1.0520338 

-2.3720509 

-1.615E-03 

-3.777E-03 

1.990E-03 

8.921E-04 

£x 


2.53E-03 


4.65E-03 

dui/dp 55.4605800 
du 3 /dp -56.5052230 

55.4057600 

-56.4808140 

9.884E-04 

4.320E-04 

55.3658650 

-56.9128710 

1.708E-03 

-7.214E-03 

£u 


1.08E-03 


7.41E-03 

dg2/dp 0.67758780 

0.6766357 

1.405E-03 

0.6767562 

1.227E-03 


eg 1.41E-03 1.23E-03 

Table 5.16 Central Difference Approximations to the Parameter Sensitivities for problem 3 
parameter 3 


Kuhn Tucker 

Forward Difference Approximations 


Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp -10.0000000 

-9.99999986 

1.4E-08 

-9.99999986 

1.4E-08 

dxi/dp 1.3057920 
d X2 /dp 0.5460588 
dx 3 /dp -1.0541310 
dx 4 dp -2.3741690 

1.3126912 

0.5453465 

-0.9886692 

-2.3442141 

-5.124E-02 

1.304E-03 

6.210E-02 

1.262E-02 

1.3151920 

0.5489030 

-1.0457860 

-2.3672333 

-7.199E-03 

-5.209E-03 

7.916E-03 

2.921E-03 

£x 


8.15E-02 


1.23E-02 

dui/dp 55.4605800 
du 3 /dp -56.5052230 

56.8153210 

-56.9128710 

-2.443E-02 

-7.214E-03 

55.5101970 

-56.7627370 

-8.946E-04 

-4.557E-03 

£u 


2.55E-02 


4.64E-03 

d g2 /dp 0.67758780 

0.6668760 

1.581E-02 

0.6757969 

2.643E-03 

e g 


1.58E-02 


2.64E-03 


Table 5.17 Forward Difference Approximations to the Parameter Sensitivities for problem 3 
parameter 3 


5.3.4 Problem 4 

The final problem in the test set contains 4 parameters, the estimated sensitivities are 
presented in Tables 5.18 - 5.23. As before, the central difference approximations produce 
the best results, and for this problem, the small step perturbation performs best. Excellent 
agreement was achieved for all the parameters in this problem. 
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Kuhn Tucker 

Central Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

0.517404073 

0.51740407 

0.00E+00 

0.51740407 

0.00E+00 

dxi/dp 

-0.40000000 

-0.4000000 

0.000E+00 

-0.400000 

0.000E+00 

dx^/dp 

0.09729967 

0.0972994 

2.960E-06 

0.097300 

-3.710E-06 

dx3/dp 

-0.20000000 

-0.2000000 

0.000E+00 

-0.200000 

0.000E+00 

dxVdp 

0.28602513 

0.2860244 

2.447E-06 

0.286026 

-3.077E-06 

dx5/dp 

-0.06773997 

-0.0677393 

9.994E-06 

-0.067741 

-1.247E-05 

e x 



1.07E-05 


1.34E-05 

du3/dp 

0.40709971 

0.4071024 

-6.681E-06 

0.407102 

-6.731E-06 

dus/dp 

-0.05662417 

-0.0566248 

-1.038E-05 

-0.056625 

-1.038E-05 

du6/dp 

0.34730309 

0.3473052 

-6.191E-06 

0.347305 

-6.219E-06 

duVdp 

0.01908492 

0.0190855 

-3.280E-05 

0.019086 

-3.233E-05 

£u 



3.56E-05 


3.52E-05 


Table 5.18 Central Difference Approximations to the Parameter Sensitivities for problem 4 
parameter 1 



Kuhn Tucker 

Forward Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

0.517404073 

0.51740407 

0.00E+00 

0.51740407 

0.00E+00 

dxi/dp 

-0.40000000 

-0.4000000 

O.OOOE+OO 

-0.400000 

0.000E+00 

dx2/dp 

0.09729967 

0.0973210 

-2.193E-04 

0.097289 

1.053E-04 

dx3/dp 

-0.20000000 

-0.2000000 

0.000E+00 

-0.200000 

0.000E+00 

dx4/dp 

0.28602513 

0.2860770 

-1.814E-04 

0.286000 

8.713E-05 

dxs/dp 

-0.06773997 

-0.0677900 

-7.388E-04 

-0.067716 

3.548E-04 

e x 



7.92E-04 


3.80E-04 

du 3 /dp 

0.40709971 

0.4073857 

-7.025E-04 

0.407116 

-4.073E-05 

dus/dp 

-0.05662417 

-0.0565732 

8.997E-04 

-0.056623 

2.283E-05 

du6/dp 

0.34730309 

0.3475276 

-6.464E-04 

0.347315 

-3.481E-05 

du9/dp 

0.01908492 

0.0190319 

2.780E-03 

0.019083 

9.583E-05 

£u 



3.07E-03 


1.12E-04 


Table 5.19 Forward Difference Approximations to the Parameter Sensitivities for problem 4 
parameter 1 
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Kuhn Tucker Central Difference Approximations 


Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 0.306110869 

0.30611087 

0.00E+00 

0.30611087 

O.OOE+OO 

d Xl /dp 0.00000000 
d X2 /dp -0.14711868 
dx 3 /dp 0.00000000 
dxVdp -0.04916518 
dx 5 /dp 0.09817962 

0.0000000 

-0.1471186 

0.0000000 

-0.0491649 

0.0981794 

0.000E+00 

6.797E-07 

0.000E+00 

4.841E-06 

2.332E-06 

0.000000 

-0.147119 

0.000000 

-0.049165 

0.098179 

0.000E+00 

4.758E-07 

0.000E+00 

3.620E-06 

1.742E-06 

Ex 


5.42E-06 


4.05E-06 

du 3 /dp -0.05662417 
du 5 /dp 0.05051545 
d U6 /dp -0.06156394 
du 9 /dp 0.00240162 

-0.0566248 

0.0505156 

-0.0615644 

0.0024015 

-1.063E-05 

-3.662E-06 

-6.903E-06 

3.935E-05 

-0.056625 

0.050516 

-0.061564 

0.002402 

-1.038E-05 

-5.424E-06 

-7.115E-06 

-2.373E-06 

£u 


4.15E-05 


1.39E-05 


Table 5.20 Central Difference Approximations to the Parameter Sensitivities for problem 4 
parameter 2 


Kuhn Tucker Forward Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

0.306110869 

0.30611087 

0.00E+00 

0.30611087 

0.00E+00 

dxi/dp 

0.00000000 

0.0000000 

0.000E+00 

0.000000 

0.000E+00 

dx2/dp 

-0.14711868 

-0.1470723 

3.150E-04 

-0.147117 

1.196E-05 

dx 3 /dp 

0.00000000 

0.0000000 

O.OOOE+OO 

0.000000 

0.000E+00 

dx4/dp 

-0.04916518 

-0.0490525 

2.292E-03 

-0.049161 

8.675E-05 

dx 5 /dp 

e x 

0.09817962 

0.0980709 

1.107E-03 

2.56E-03 

0.098176 

4.190E-05 

9.71E-05 

du 3 /dp 

-0.05662417 

-0.0570971 

-8.352E-03 

-0.056648 

-4.279E-04 

dus/dp 

0.05051545 

0.0510986 

-1.154E-02 

0.050545 

-5.817E-04 

du^dp 

-0.06156394 

-0.0620483 

-7.868E-03 

-0.061589 

-4.015E-04 

du 9 /dp 

Eu 

0.00240162 

0.0024686 

-2.790E-02 

3.23E-02 

0.002405 

-1.404E-03 

1.63E-03 


Table 5.21 Forward Difference Approximations to the Parameter Sensitivities for problem 4 
parameter 2 
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Kuhn Tucker Central Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

1.183954566 

1.18395457 

0.00E+00 

1.18395457 

0.00E+00 

dxi/dp 

-0.20000000 

-0.2000000 

O.OOOE+OO 

-0.200000 

0.000E+00 

dx2/dp 

0.09408231 

0.0940821 

2.689E-06 

0.094082 

2.317E-06 

dx3/dp 

-0.35000000 

-0.3500000 

0.000E+00 

-0.350000 

0.000E+00 

dxVdp 

0.22881747 

0.2288169 

2.666E-06 

0.228817 

2.273E-06 

dx 5 /dp 

£x 

0.02931310 

0.0293137 

-2.030E-05 
2. 06 E -05 

0.029314 

-1.723E-05 

1.75E-05 

du3/dp 

0.34730308 

0.3473054 

-6.709E-06 

0.347305 

-6.277E-06 

dus/dp 

-0.06156394 

-0.0615644 

-8.089E-06 

-0.061564 

-7.115E-06 

dus/dp 

0.72069515 

0.7206968 

-2.248E-06 

0.720697 

-2.456E-06 

dus/dp 

Eu 

0.15964688 

0.1596474 

-3.382E-06 

1.13E-05 

0.159647 

-3.570E-06 

1.04E-05 


Table 5.22 Central Difference Approximations to the Parameter Sensitivities for problem 4 
parameter 3 


Kuhn Tucker Forward Difference Approximations 


Method 

AP =2% 

relative error 

AP=0.1% 

relative error 

df/dp 1.183954566 

1.183954566 

0.00E+00 

1.18395457 

0.00E+00 

dxj/dp -0.20000000 
d X2 /dp 0.09408231 
dx 3 /dp -0.35000000 
dxVdp 0.22881747 
dx 5 /dp 0.02931310 

-0.2000000 

0.0941409 

-0.3500000 

0.2289599 

0.0291757 

0.000E+00 

-6.226E-04 

0.000E+00 

-6.226E-04 

4.687E-03 

-0.200000 

0.094083 

-0.350000 

0.228818 

0.029313 

0.000E+00 

-2.073E-06 

0.000E+00 

-2.098E-06 

1.559E-05 

e x 


4.77E-03 


1.59E-05 

du 3 /dp 0.34730308 
du 5 /dp -0.06156394 
du 5 /dp 0.72069515 
du 5 /dp 0.15964688 

0.3485794 

-0.0614102 

0.7229895 

0.1595192 

-3.675E-03 

2.497E-03 

-3.184E-03 

8.000E-04 

0.347369 

-0.061557 

0.720811 

0.159641 

-1.895E-04 

1.153E-04 

-1.611E-04 

3.620E-05 

Eu 


5.52E-03 


2.77E-04 


Table 5.23 Forward Difference Approximations to the Parameter Sensitivities for problem 4 
parameter 3 
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Kuhn Tucker 

Central Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

0.010389619 

0.01038962 

0.00E+00 

0.01038962 

0.00E+00 

dxi/dp 

dx2/dp 

dx3/dp 

dxVdp 

dxs/dp 

0.00000000 

-0.03975934 

0.00000000 

0.07614087 

0.15499104 

0.0000000 

-0.0397594 

0.0000000 

0.0761408 

0.1549911 

0.000E+00 

-8.803E-07 

0.000E+00 

1.116E-06 

-5.162E-07 

0.000000 

-0.039759 

0.000000 

0.076141 

0.154991 

0.000E+00 

-1.082E-06 

0.000E+00 

1.366E-06 

-6.452E-07 

e x 



1.51E-06 


1.86E-06 

du3/dp 

dus/dp 

du6/dp 

du9/dp 

0.01908492 

0.00240162 

0.15964687 

0.08386033 

0.0190855 

0.0024016 

0.1596474 

0.0838607 

-2.950E-05 

-7.745E-06 

-3.382E-06 

-4.031E-06 

0.019086 

0.002402 

0.159647 

0.083861 

-3.238E-05 

-2.498E-06 

-3.633E-06 

-4.042E-06 

6u 



3.09E-05 


3.29E-05 


Table 5.24 Central Difference Approximations to the Parameter Sensitivities for problem 4 
parameter 4 



Kuhn Tucker 

Forward Difference Approximations 



Method 

AP =2% 

relative error 

AP=0. 1 % 

relative error 

df/dp 

0.010389619 

0.01038962 

0.00E+00 

0.01038962 

0.00E+00 

dxi/dp 

dx^dp 

dx3/dp 

dx^dp 

dxs/dp 

0.00000000 

-0.03975934 

0.00000000 

0.07614087 

0.15499104 

0.0000000 

-0.0397539 

0.0000000 

0.0761542 

0.1549782 

0.000E+00 

1.374E-04 

0.000E+00 

-1.744E-04 

8.265E-05 

0.000000 

-0.039759 

0.000000 

0.076143 

0.154989 

0.000E+00 

1.856E-05 

0.000E+00 

-2.359E-05 

1.123E-05 

£x 



2.37E-04 


3.20E-05 

du3/dp 

dus/dp 

du^/dp 

du9/dp 

0.01908492 

0.00240162 

0.15964687 

0.08386033 

0.0193233 

0.0024360 

0.1599642 

0.0841484 

-1.249E-02 

-1.432E-02 

-1.988E-03 

-3.435E-03 

0.019097 

0.002403 

0.159663 

0.083875 

-6.561E-04 

-7.328E-04 

-1.032E-04 

-1.755E-04 

£u 



1.94E-02 


1.00E-03 


Table 5.25 Forward Difference Approximations to the Parameter Sensitivities for problem 4 parameter 4 


5.3.5 A Graphical Comparison of the Accuracy Delivered by Several Different 
Methods 

This section presents a graphical method that can be used to assess the accuracy of 
the parameter sensitivity derivatives calculated by various methods. When we plot a bar 
chart of the calculated sensitivity derivatives versus the true sensitivity derivatives, we can 
identify which components of the gradient are least accurate. A brief discussion of the 
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significance of each graph will be provided. 


Figure 5.1 and 5.2 present a graphical comparison of the dx/dpi and 3u/3pi (as 
calculated by various methods 2 ) for problem 3. In addition to results discussed in section 
5.3.3, we calculated sensitivities by the modified Kuhn-Tucker method. The Hessian 
approximations are provided in section 5.2. 


We can see in figure 5.1 that for large values of Ap for forward differencing we did 
not obtain results as accurate as those delivered by the other variants of the RQP sensitivity 
algorithm. If we examine the bars for the KT w/ BFS Hess we see that there is some 
discrepancy in the calculated values of the sensitivity derivatives. We also can see that 
using the Hessian approximation delivered by using the SRI update we were able to obtain 
results better than those obtained by using the BFS approximation. 



■ KT 

^ RQP cd AP=2% 
H RQP cd AP=0.1% 
m RQP fd AP=2% 
□ RQPfd AP=0.1% 
CD KTw/ BFS Hess 
U KT w/SRl Hess 


Figure 5.1 A comparison of the accuracy of dx/dp for Problem 3 parameter 1 


We can see in figure 5.2 that all but the KT w/ BFS Hess method produces good 
estimates of the sensitivity derivatives. The discrepancy in 3u/3p is due to the Hessian 

2 KT - Kuhn Tucker sensitivity equations with exact derivatives 

RQP cd AP=2% - RQP algorithm using 2 iterations to solve perturbed problem with central difference 
approximations Ap=2% of the nominal value 

RQP cd AP=0.1% - RQP algorithm using 2 iterations to solve perturbed problem with central difference 
approximations Ap=0.1% of the nominal value 

RQP fd AP=2% - RQP algorithm using 2 iterations to solve perturbed problem with forward difference 
approximations Ap=2% of the nominal value 

RQP fd AP=0.1% - RQP algorithm using 2 iterations to solve perturbed problem with forward difference 
approximations Ap=0.1% of the nominal value 

KT w/ BFS Hess - The Kuhn Tucker sensitivity equations were solved using the approximate Hessian 
delivered by RQOPT when using the BFS update 

KT w/ SRI Hess - The Kuhn Tucker sensitivity equations were solved using the approximate Hessian 
delivered by RQOPT when using the SRI update 


64 





approximation not converging. 
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Figure 5.2 A comparison of the accuracy of 8u/9p for problem 3 parameter 1 


Figures 5.3 and 5.4 present comparisons of the accuracy of the sensitivity 
derivatives (calculated by various methods 2 ). We can see that all methods are in good 
agreement with the exception of the calculation of the KT with the approximate Hessian. 
Again the Hessian approximation did not converge fully for this problem. 
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Figure 5.3 A comparison of dx/8p for Problem 4 as calculated by various methods 
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Figure 5.4 A comparison of various methods for calculation of du/dp for Problem 4 


5.4. CONCLUSIONS FROM TEST RESULTS 


In summary, we feel that this initial testing has shown that the RQP based method 
can produce reliable estimates of parameter sensitivities. Further testing is required to 
determine its performance on engineering problems, and to further resolve questions about 
algorithm parameters and variations. The following conclusions can also be drawn from 
this initial testing. 

We saw in section 5.2 that the Hessian approximation improved if we allowed the 
approximation to be updated during the sensitivity analysis. During testing of the other 
problems in the test set we observed the Hessian approximation improving as we calculated 
the parameter sensitivities. This implies that if the Hessian approximation did not converge 
during the solution of the original problem (or converged in a projected sense) then a good 
estimate of the Hessian approximation can be built during the sensitivity analysis. Once we 
have a good approximation of the Hessian approximation then we can switch from a more 
expensive central differencing approximation to a less expensive forward difference 
approximation to the sensitivity derivative approximation. This conclusion is also 
encouraging because the RQP sensitivity algorithm approaches the Kuhn-Tucker sensitivity 
algorithm as the approximation improve (as was demonstrated in section 3.2). 

If the Hessian approximation found by the RQP method is a good approximation of 
the Hessian of the Lagrangian then the Kuhn-Tucker sensitivity equations can be used with 
the approximate Hessian. However, tests need to be developed to check if the Hessian 
approximation has converged. If the Hessian approximation did not converge then 
updating it while using the RQP sensitivity algorithm can cause the Hessian approximation 
to converge. When this happens the Kuhn-Tucker sensitivity equations may be used with 

the Hessian approximation. Again, the issue to be resolved before this option can be 
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investigated is how to test a Hessian approximation for convergence. 

Our experiments with the SRI update were very encouraging. Section 5.2 
demonstrated how the SRI update was able to deliver a more accurate estimate of the 
Hessian approximation than the BFS update. We also observed exact convergence of the 
Hessian approximation with the SRI update if we allowed the Hessian approximation to be 
updated during the sensitivity analysis. Near convergence of the Hessian approximation 
for problem 1 and 3 was also obtained when the SRI update was used. 

Even though our experiments with the SRI update were encouraging the SRI 
update has some serious drawbacks that will have to be investigated further before the 
method can be recommended. Unlike the BFS update the SRI update is not self-correcting 
(Ip 1987). As we saw in problem 4, the SRI update delivered an Hessian approximation 
that was singular, and the RQP method performs best if the Hessian approximation is 
positive definite. 

In section 5.3, we observed that the method was able to approximate sensitivity 
derivatives. We saw that using central differencing yielded more accurate estimates of the 
sensitivity derivatives than using forward differencing approximations. We did not 
observe any major sensitivity to Ap when we were using the central differencing option. 

We also observed that the forward differencing approximation has more trouble 
evaluating sensitivity derivatives, and the perturbation step size Ap had more effect when 
the functions being approximated are nonlinear in the parameter. 

Finally, we observed that when the sensitivity derivative is small, the RQP 
sensitivity algorithm can sometimes have a difficult time finding the exact value. This was 
demonstrated in the calculation of dnjjdp for problem 1 and also in the calculation of the 
sensitivity derivatives for p 2 in problem 3. 
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6 . 


Changes in the Active Set 


This chapter describes our findings regarding changes in the active set of 
constraints. The first section describes the four cases of changes in the active set of 
constraints. The second section provides sample plots illustrating the predicted behavior 
from section 1. Next, the third section discusses several proposed solutions for dealing 
with changes in the active set. The final section provides some procedures that can be used 
to generate problems where there are changes in the active set. 


6. 1 . CASES OF ACTIVE SET CHANGE 

When the active set of constraints change some of the sensitivity derivatives may be 
discontinuous. It is also possible that at the point where the active set changes, the 
gradients of the constraints may become linearly dependent. When this happens the 
optimization problem may become very difficult to solve. The four cases for changes in the 
active set are; 

1 . A constraint enters the active set (constraint gradients are linearly independent). 

2. A constraint leaves the active set. 

3. A constraint enters the active set (constraint gradients are linearly dependent). 

4. A constraint enters the active set and causes there to be no feasible solution 

Case 1 and Case 2 are complementary 1 . This can best be demonstrated by an 
example. Consider an optimization problem whose solution changes as a parameter p, is 
varied and the active set changes as p increases. Assume the original problem is optimized 
for a p=pO and and a sensitivity analysis is performed. The value of p is then increased to 
p=pl (where pi > pO) and the problem is reoptimized. At p=pl, constraint g 3 will enter the 
active set if p is increased any further, but the gradients of the constraints remain linearly 
independent as p is increased. Now p is increased to p=p 2 (where p 2 > pi > pO) and at this 
point constraint g 3 has been added to the active set. A Case 1 active set change has 
occurred. To see the complementarity, assume we begin with p=p 2 and conduct a 
sensitivity analysis for decreasing values of p. Now when p = pi, for any value of p less 
than pi constraint g 3 will leave the active set. Thus, when p reaches pO, constraint g 3 will 
have left the active set and we have had a Case 2 change in the active set for p going from 
p 2 to pO. 

The behavior of the sensitivity derivatives for Cases 1 and 2 is characterized as 

1 This means that active set change algorithm that are developed for case 1 changes in the active set may 
also be used for case 2 changes in the active set 
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follows. There is a discontinuity in the rate of change of the objective function with respect 
to p (d2f*/dp2 is discontinuous), there can be a discontinuity in the rate of change of the 
some of the design variables (dx*/5p is discontinuous), and the rate of change of the 
Lagrange multipliers can be discontinuous (i.e. 3u*/dp for a constraint leaving the active 
set will go from some nonzero value to zero, this may cause a change in the rate of change 
of the other Lagrange multipliers as well). The Hessian of the Lagrangian is continuous 
with respect to variations in p. For these problems there will be directional derivatives of 
dx*/dp and du*/dp that can be used to make estimates of how the design variables and 
Lagrange multipliers will change. Second order extrapolations of the behavior of the 
objective function can also be made using d2f*/dp2 in a directional sense. 

The characteristic of Case 3 is a discontinuity in the Lagrange multiplier estimate 
(u* is discontinuous). The discontinuity in the Lagrange multiplier causes a discontinuity 
in df*/dp and also causes the Hessian of the Lagrangian to be discontinuous. Since the 
active set changes, there will also be a discontinuity in (dx*/9p). At the point where the 
constraints become linearly dependent it will not be possible to solve the Kuhn-Tucker 
sensitivity equations because they will become singular. For Case 3, there may be an 
exchange of constraints in the active set (i.e. the new constraint may replace one of the 
constraints that is already in the active set as p moves through the point). 

The main characteristic of Case 4 is that there only exists a directional derivative 
away from the point where the path terminates. In order to calculate the directional 
derivative for Case 4 we need to be able to find the proper active set to follow when we 
leave the degenerate point. Case 4 changes in the active set will be common when the user 
overconstrains the design, i.e. sets the performance specifications to high to be physically 
meet by the given design configuration. 


6.2. DEMONSTRATION OF CASES 

This section describes the behavior that we observed in the initial test set for 
problems that had changes in the active set. This section also discusses two test problems 
that are not in the initial test set that are used to demonstrate the behavior of Case 3 changes 
in the active set. 

The problems in the initial test set only had Case 1 and Case 2 changes in the active 
set. Indications of where the active set changes were shown on the plots of the optimum 
sensitivity provided in the appendix. 

If we examine the plots that are presented in the appendix (figures A.1-A.8, A.10, 
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and A.l 1) we can see the discontinuity in 9x*/9p and 9u*/9p when the active set changes. 
In problems 3 and 4, because there is only a small discontinuity in dx/dp, it is difficult to 
see in the graphs provided (figures A.5,6,7,10,1 1) in the appendix. However for 
problems 3 and 4 it is easy to see the discontinuity in 9u/9p, the slope of the Lagrange 
multipliers, when the active set of constraints changes. 

If we examine figures A.4 and A. 11, we see that not all of the components of 
9x*/9p are discontinuous when the active set changes. For problem 2 in figure A.4 we do 
not see a discontinuity in 9x2*/9p, and for problem 4 in figure A. 10 we do not see a 
discontinuity in 9xi*/9p or 9x3*/9p, when the active set changes. This behavior is partially 
due to the fact that we were studying right hand side perturbations for these problems. The 
symmetry in problem 2 can be used to explain its behavior. In problem 4, the gradient of 
constraint g 9 had zero’s in the first and third locations, thus we expect that perturbations of 
constraint g 9 should have no effect on 9xi*/9p or 9x3*/9p. 

The discontinuity in d 2 f*/dp 2 is difficult to see, when the active set changes for 
most of the problems in the initial test set. However for test problem 2, parameter 1, 

(figure A.3) we can see the discontinuity in d 2 f*/dp 2 . In the plot of the optimum value of 
the objective function versus pi we can see that for values of pi < 1.055, d 2 f*/dp 2 is less 
than zero (a region of negative curvature) and for values of pi > 1.055, d 2 f*/dp 2 is greater 
than zero (a region of positive curvature). 

In problem 2 (p 2 ), problem 3 (pi,p2,P3) and problem 4 (p 4 ) we studied the effect of 
perturbing the right hand sides of the inequality constraints. In these problems we studied 
a range of perturbations from where the constraint was inactive to where the constraint is in 
the active set. When the constraint was inactive perturbations in the parameters had no 
effect on the optimum. When the constraints entered the active set for these problems they 
caused d 2 f*/dp 2 to go from zero to some nonzero number. The discontinuity of d 2 f*/dp 2 
can be seen in figures A.4, A. 5, A.6, A.7, and A.l 1. 

The following example will demonstrate how the optimum behaves for Case 3, 
when the gradients of the constraints become linearly dependent upon a change in the active 
set of constraints. 

minimize: f = xi2 + (P-l)2 (6.1) 

subject to: gi = 3xi +2 P- 10 >0 (6.2) 

g2 = 2 xi + 3 P - 10 > 0 (6.3) 

when P = 2, the minimum f* = 5 occurs at xi* = 2 with both constraints active and the 

gradients of the constraints linearly dependent. The Lagrange multipliers and in the family 
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ui>U2 e {3 ui + 2 U 2 = 4, ui > 0, U 2 > 0} (6.4) 

For this problem df*/dp, 5x*/9p and du/9p do not exist. These derivatives do exist in 
both the positive and negative directions and are indicated by dx/dp + for increasing values 
of p and dx/dp- for decreasing values of p. 

Figure 6.1 presents the sensitivity information for problem 6. 1-6.3. Figure 6.1 (a) 
and (b) represents the first order predictions of the new values of the Lagrange multipliers 
for this problem. For this problem the linear predictions agree with the optimum Lagrange 
multipliers. There is a discontinuity at Ap = 0.0, therefore there are only be directional 
derivatives for these values. Figure 6.1 (c) represents linear predictions of the new value 
of the objective function. Notice again that there is a discontinuity in the slope of the 
prediction. Thus df*/dp does not exist for this value of p (however df*/dp + and df*/dp* 
will exist in a directional sense where the superscript + or - indicate the direction of change 
in p). Figure 6.1 (d) represents the predicted location of the optimum. Notice the 
discontinuity in the slope at Ap = 0.0, the true location of the optimum agrees with the 
linear prediction for this problem. 



Figure 6.1 Plot of the optimum sensitivity predictions for problem 6.27-29 

We have also examined a test problem used by Fiacco and Ghaemi (1982). This 
test problem has a linear dependence in the constraints gradients. The problem is to design 
a corrugated bulkhead for an oil tanker. The objective function is to minimize the weight of 
the bulkhead with constraints placed on allowable bending stress, moment of inertia, 
corrosion, and minimum gage thickness of the material. The problem has six design 
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variables, 18 design parameters, 17 constraints, a nonlinear objective function, and near 
linear constraints. We did not apply any scaling to this problem, even though the 
constraints for this problem are poorly scaled. 

A sensitivity analysis was performed for variations in the parameter, LT, from 
0.01% to 1% using the RQP code RQOPT (Beltracchi and Gabriele 1987). One step 
convergence of the objective function value, predicted objective function value, design 
variables, and the Lagrange multipliers were also used. This sensitivity analysis was 
performed before the RQSEN system was built, and all of the calculations were performed 
by hand. Therefore we did not record all of the values of the design variables, Lagrange 
multipliers, and constraints. 

Figure 6.2(a) is a plot of the sensitivity of the optimal objective function for the 
bulkhead problem versus the parameter LT which was the length to the top brace. At LT = 
482.8 we notice a discontinuity in the slope of the objective function. At this point we have 
encountered a Case 3 change in the active set. At LT = 476.25 there is a Case 2 change in 
the active set. 

Figure 6.4 presents a plot of the value of constraint 13 (g^) versus the value of LT. 
At LT~ 482.8, the value of this constraint goes to zero and the constraint enters the active 
set. 


Figure 6.3 presents a plot of the value of the Lagrange multiplier for constraint 12 
versus LT. Notice that there is a discontinuity at LT = 482.8, this is due to a linear 
dependence of the constraint gradients of the active set. As LT moves from one side of 
482.8 to the other, there will be a change in the active set of constraints with gj 3 replacing a 
constraint on the lower bound of xg. 

Figure 6.2(d) presents a plot of the optimal value of design variable six. For LT > 
482.8, the lower bound constraint is active. However, when LT < 482.8 the optimal value 
of this design variable is no longer on the bound. Figure 6.2 also presents plots of the 
behavior of xj and X 3 , we can also see discontinuities in these variables when the active set 
changes. 
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Figure 6.2 Sensitivity of the Bulkhead problem 
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Figure 6.3 Plot of The Lagrange Multipliers for the Bulkhead problem 



Figure 6.4 Plot of gi 3 versus Lt for the Bulkhead problem 


For this problem there is a second change in the active set at LT ~ 476.25 when gj 2 
leaves the active set. As LT moves through 476.25 there is not the same type of 
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discontinuity as we saw for Case 3. This represents Case 1, a constraint leaving the active 
set. We can see some of the characteristics of this case in the plots. In figure 6.2, the 
slope of the objective function is continuous at LT = 476.25. In figure 6.3, the value of the 
Lagrange multiplier of constraint g^ goes to zero signifying the constraint has left the 
active set. Finally, in figure 6.2 we can see that dxj/dp, dxydp, and dx^/dp are 
discontinuous at LT = 476.25. 

A Case 4 change in the active set can be illustrated by the following simple one 
variable problem. 

minimize f(x) = x 

Subject to: gi = x - pi > 0 
g2 = x - 2 > 0 

When pi = 1 the optimum solution to this problem is f*=1.0, x*=1.0, and constraint gi is 
active. When pi is increased to pi=2.0 constraint g2 enters the active set. However, if pi 
is increased beyond 2.0 there is no feasible solution for this problem. 

An engineering example of a case 4 change in the active set can be illustrated by the 
following example. Find the optimum air plane to fly a given mission, where the design 
variables may be the size of the wings, engine size, cruise altitude, etc. The design 
parameters might be; total cargo weight, runway length, air temperature at take-off, etc. 
Assume a constraint on take-off distance is applied, 

Take-off distance < Available Runway Length 
and a parameter sensitivity study of total cargo weight is performed. As the total cargo 
weight is increased, the take-off distance of the airplane increases, and the design variables 
of the airplane may also change. If the total cargo weight is increased beyond a certain limit 
it may become impossible for the aircraft to take-off at the given runway length. A 
parameter sensitivity analysis for values of the total cargo weight for values greater than 
this limit are meaningless because there will be no solution to the problem since the plane 
cannot take-off. Thus, for this problem when the runway length constraint enters the active 
set there will be a case 4 change in the active set causing no feasible solution to exist. 


6.3. PROPOSED SOLUTIONS 

This section will discuss some algorithms that can be used to obtain more accurate 
estimates of the parameter sensitivity after the active set of constraints has changed. The 
first subsection describes algorithms that can be used to calculate sensitivity derivatives 
when there are Case 1 or Case 2 change in the active set. The second subsection presents 
an example demonstrating how the algorithms presented in the previous section work. The 
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final subsection discusses how sensitivity derivatives can be calculated for Case 3 changes 
in the active set. 

As we have mentioned previously, the sensitivity derivatives do not exist at points 
where the active set changes due to the discontinuities present in the optimum path. 
However, the sensitivity derivatives will exist in a directional sense. The proposed 
algorithms for dealing with changes in the active set are based on calculating directional 
derivatives. 

Directional derivatives can be approximated as follows 


df*(Ap + ) Limit f f*(pO + Ap) - f*(pO)> 

~dp ~Ap— >0 + |^ ^ 


(6.5) 


Here df*(Ap + )/dp indicates the rate of change of a function when a parameter is perturbed 
in the positive direction. When the base point for a sensitivity analysis is degenerate, we 
can approximate directional derivatives using the RQP sensitivity algorithm in the same 
way as we approximate other derivatives. 


6.3.1 Dealing with Case 1 and 2 


We begin this section by presenting the following example problem 


minimize: f = (xj + 1)2 + (x 2 - 2)2 (6.6) 

subject to: gi = xi - p > 0 (6.7) 

g 2 = 6-2xi-x 2 >0 (6.8) 


When p = pO = 1 as shown in figure 6.5, f* = 4, x* = (1,2). A sensitivity analysis 
indicates that as p is increased, 8x*/8p = (1,0). Notice, that as we increase p, eventually 
constraint g 2 will become active which is a Case 1 change in the active set. Increasing p 
from (x*,pl) will result in the optimum point moving along the intersection of the two 
constraints. The deflection algorithm is based on finding a constraint intersection and then 
calculating new values of 9x*/9p and 9u*/9p along the new active set. 
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Figure 6.5 Deflection of the Sensitivity Search Direction 

The proposed algorithm for dealing with constraints entering or leaving the active 
set is given by the following steps: 


1 . Determine 3x*/9p at the base point 

2. Calculate Apl for intersection of the constraint by finding the minimum Ap that 
causes a constraint to enter/leave the active set 


mm 


Si 


min ( u; \ . 

- 1 — j € Active set 

11. X J 


j € Active set 



(6.9) 

( 6 . 10 ) 


3. 

4. 


Calculate dx* V9p with the active set updated to reflect the change indicated by 
step 2. 

Calculate Ax and Au by 


dx* 

a * = %t a p 


jp 

a du* * 

A U =^-Ap 


1 + 
1 + 


ax*i 

V Ap2 

3u*l A O 

"35" A P 2 


( 6 . 11 ) 

( 6 . 12 ) 


where 
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Apl = min ( p - pO, 0) 
Ap2 = max ( p - pi, 0) 


(6.13) 

(6.14) 


5. Calculate the new value of the objective function by 

fnew = f(x*,p0) + Ap ^p-+|^-Ax (6.15) 

The first step in the above algorithm is to perform a sensitivity analysis to determine 
which direction the optimum will move in. In the second step, we attempt to identify the 
value of p for which the active set will change using equations 6.9 and 6.10. A new search 
direction is calculated in step three which includes the changed active set. In step four, a 
formula is provided which defines the proper value for Ax for before the constraint is 
encountered and after the constraint becomes active (inactive). This result is used in step 5 
to predict a new value for the optimum at the new point. 

In step three of the deflection algorithm, it is necessary to find the deflected search 
direction at the point where the new constraint enters the active set. This can be done by 
adding the newly violated constraint (or removing the constraint that is leaving the active 
set) to the Kuhn-Tucker sensitivity equations and then solving for dxV5p and 9uVdp. 

This computation can be performed efficiently by using matrix updating techniques as 
described in (Diewart 1984). In order to obtain a more accurate estimate of the 
sensitivities, the gradients of the active constraints can be re-evaluated at x = xl, where x 1 
is the predicted intersection of the new constraint. 

When using the deflection algorithm, a check should be made to assure that the 
constraint gradients are not linearly dependent. If the constraint gradients are linearly 
dependent then this procedure cannot be used to solve for the optimum sensitivities because 
there is no solution to the Kuhn-Tucker sensitivity equations. The procedure to be used 
when the constraint gradients are linearly dependent will be discussed in section 6.3.3. 

The estimate obtained in step 5 is a linear extrapolation. A better estimate of the 
new value of the optimum can be made using a quadratic extrapolation. The formula for a 
quadratic extrapolation is (McKeown 1980, Fiacco 1983, Barthelemy and Sobieski 1983) 

df* 1 d2f* 

f*new = f(x*) + Api J Api Api (6. 1 6) 

where 

d 2 f* a 2 L r a 2 L it9x* auTag avTah 

^T’35T 3£ + 3FT 35? (6 ' 17) 
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Using (6.16) to predict the new value of the objective function as p varies yields a 
more accurate estimate of the new optimal objective function until the active set changes. 
We should note that equation 6.17 is in terms of dg/dp and dh/dp, and can be rewritten in 
the following form in terms of dx*/dp 


4 , 92l ~ | t dx* dx*T d 2 L dx* 
dpi 2 dpj 2 + Upi^x*] ^pT + ^Pi" dx* 2 "5pT 


(6.18) 


If we substitute equation 6.18 into equation 6.16 and write the equation in terms of 
Ap and Ax we will obtain 

« * 1T A d 2 L A „ a d 2 L T a a t d 2 L A 1 

f*new = fisto rder + j|^Api ^2 A Pi + 2Api^p-Ax + Ax T ^-AxJ (6.19) 

where the first order approximation comes from equation 6.15. This form has the 
advantage of being able to predict the new value of the objective function at a new value of 
x and p. When the deflection algorithm is used to calculate the new location of the 
optimum then equation 6.19 can be used to make a second order estimate of the new value 
of the objective function. 


As mentioned in section 6.1, a constraint leaving the active set (Case 2) is 
complementary to a constraint entering the active set (Case 1). It may be possible to predict 
df*/dp when a constraint leaves the active set without calculating dx ! */dp and duVdp along 
the new active set. This could be beneficial for instances that only require estimates of 
f*(p). To predict the behavior of the optimum when a constraint leaves the active set we 
can use directional derivatives. Using the following formula (Fiacco 1983) 


df* df dhi 

dF = 3p + Ji Vi 3p 




( 6 . 20 ) 


and the fact when an inequality constraint leaves the active set its Lagrange multiplier goes 
to zero, we can obtain an estimate of df*/dp at the point where the change in the active set 
takes place, this is done by obtaining an estimate of du/dp from a sensitivity analysis at the 
base point using it be used to estimate when a constraint will leave the active set. A 
prediction of the df*/dp when the active set changes can then be made from the following 
formula 


df* df ldhj 



where 



( 6 . 21 ) 


( 6 . 22 ) 


79 



(6.23) 


u i = Ui * + 3jr Apl 

When these values of u 1 and v 1 are used we obtain a new estimate for the value of the 
df*/dp that will be valid when the value of p is increased or decreased past p 1 . 

The deflection algorithm can be used to predict when a second change in the active 
set will take place. A prediction of when a second constraint may be added to or removed 
from the active set can be made by making a linear approximation using the following 
formulas 



dx 


A 1 ' 

APi J 


g Active set 


(6.24) 


Uj =Uj 



e Active set 


(6.25) 


and 


Apf = 


min 


VvPi OX* dpi J 


j g Active set 


Apj =min 


f “j A 


r*ai' 

vj_ dpi . 


J 


j e Active set 


(6.26) 


(6.27) 


where Ap 2 predicts the value of pi that will cause the second constraint to enter (using eq. 

6.26) or leave (using eq. 6.27) the active set, gj is a prediction of the value of the j* 

constraint, and Uj 1 is the predicted value of the Lagrange multiplier when the active set 
changes. The predicted value for Ap2 is calculated as a linear approximation of the value 
of the constraint in the dxVdp direction. If the constraints are interrelated (i.e. are 
evaluated as a set) then it may be possible to use a more accurate estimate of dg/dx in 
formula 6.26, by predicting its new value by using the formula 

+ (6.28) 

When the constraints are interrelated d 2 g/dx*dp can be evaluated when dV x L/dp is 
evaluated. 
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If the calculation of the gradient of a particular constraint is not expensive then the 
value of the new gradient may be used in formula 6.26. 

The deflection algorithm and the variants that were introduced in this section will be 
effected by any nonlinearity that is present in the problem, particularly nonlinearity in the 
constraints. It should be emphasized that these sensitivities are only estimates of how the 
sensitivity will change. 

6.3.2 An Example 

This section presents an example problem to demonstrate how the deflection 
algorithm of section 6.2.1 performs. Example problem (6.6-8) will be used. Figure 6.5 
shows the solution for p = 1.0, this example assumes that p is increased to p = 2.5. The 
exact value of the new optimum at p = 2.5 is f* = 13.25, X* = (2.5,1). 

For p = p0 = 1.0, the initial search direction for an increasing p is dx*/dp = (-1,0) 
as shown in figure 6.5. Step 2 of our algorithm determines that constraint g 2 enters the 
active set when p=pl=2.0. For values of p greater than pi the location of the optimum is 
along the intersection of constraints gi and g2- The new search direction along constraints 
gl and g 2 is determined from step three to be, 

^■=(-1,2) (6.29) 

For p = p2 = 2.5, which is greater than pi, the estimated Ax is composed of the 
sum of two vectors: one vector from (x*,p0) to (x*,pl) plus the vector from (x*,pl) to 
(x*,p 2 ). Thus by equation 6.1 1 for a Ap=1.5 

Ax = (1.5,1) (6.30) 

Thus the estimate of the new location of the optimum for p=2.5 becomes 

x *new est ~ (2.5,1) (6.31) 

Which is the true location of the new optimum for this problem. Without using the 
deflection algorithm we obtain the following estimate (by equation 2.9) of the optimum 

x *new est = (2.5,2) (6.32) 


By equation 2.7 and 2.8, the new value of the objective function will be 

fKTl st = 10 (6.33) 

Using equation 6.10, which takes into account the change in the active set, we get 
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that the confidence intervals of the power plots, in particular Figure 43(c), were 
derived using the approximations of (3.30)-(3.32). Evidently, in this instance, those 
approximations combine with the nonlinear relationship between weight magnitude 
and output noise power in (4.1) to yield inadequate estimates of expected value 
and variance. 

4.3 Summary 

The expressions for expected value and variance of the array weights, output 
powers, and output power ratios have been compared with experimental results 
from a real modified SMI array which “receives” bench generated, narrowband sig- 
nals. The experimental output desired signal power, output interference power, and 
array weights all agree with the statistical predictions. The experimental output 
noise power was a few dB higher than predicted for F = 0.9. The cause of this 
problem was found to be the high variance of the second auxiliary weight which in 
turn (it was argued) was due to the relatively small input INR to that channel. In 
practice, the auxiliary elements will have approximately uniform gain across the 
field of view and thus this situation should not occur. If this problem does occur, 
it may be solved by excluding the noise eigenvector from the sample covariance 
matrix inverse as done in Section 3.3.5. 
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.4. GENERATION OF TEST PROBLEMS 


This section will describe a procedure that can be used to generate test problems that 
lossess changes in the active set. The first section describes the generation of test 
iroblems with Case 1 and Case 2 changes in the active set. The second section describes 
he generation of test problems with Case 3 changes in the active set. 

6.4. 1 Test Problems for Case 1 and Case 2 Changes in the Active Set 

Several test problems that exist in the literature possess Case 1 and 2 changes in the 
ctive set (See Schmit and Chang 1984, Vanderplaats and Yoshida 1985, Vanderplaats and 
2ai 1987). To generate new problems, active set changes for these two cases can be 
introduced into a test problem by adding constraints that are not active at the optimum but 
are violated for small changes in the parameters. The following is an outline of the steps 
that can be used to generate test problems where there are changes in the active set 

Step 1 Generate a NLP Test Problem. One such method would be to use the 
Rosen and Suzuki (1965) procedure. 

Step 2 Vary some problem parameters and find the path of the optimal solution 
with respect to p, say x*(p) from p° to p 2 . 

Step 3 Choose a value of p 1 between pO and p 2 as the point where a constraint will 
enter the active set. 

Step 4. Construct a new constraint such that 

gnew( x *(p 1 )>P 1 ) = 0 
gnew(x*(p 1 - ApXp 1 - Ap) > 0 

and the gradient of gnew(x*(p 1 ),p 1 ) is linearly independent of the gradients 
of the constraints that are already in the active set. 

Step 5. Calculate the path of the new problem from p 1 to p 2 

The above algorithm can be used to create test problems for Case 1 when a 
sensitivity analysis is conducted at p = p° and then p in perturbed to p 2 . The same test 
problem can be used as a Case 2 test problem, when the sensitivity analysis is performed at 
p = p 2 and p is moved to p°. \ 

This procedure is illustrated in figure 6.6 which shows a graph of a two variable 
test problem along with the optimum x*(pO). At the optimum, constraint gi is active. As 
the p increases from pO to pi the location of the optimum moves from x*(pO) to x*(p2), 
and constraint gi remains active. To introduce a change in the active set we can place an 
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additional constraint g 2 (shown in figure 6.6 (b)), that intersects the vector {9x(pO)/8p} Ap, 
where Ap = p2 - pO. The value of p where the path to the optimum intersects the new 
constraint will be denoted as pi. By adding constraint g 2 , the optimum location x*(p2) 
shown in figure 6.6 (b) will be different than without the constraint. 




Figure 6.6 The creation of Test Problems with changes in the Active Set 


6.4.2 Generation of Test Problems for Case 3 

This section will describe the generation of test problems that have a Case 3 change 
in the active set. We first discuss problems that are in the literature. Then we will discuss 
two different algorithms that can be used to generate test problems of this type. 

A survey of the literature revealed several test problems where the constraint 
gradients are linearly dependent when the active set changes (Case 3). Three such 
problems were found in articles by Bartholomew-Biggs (1986), Vanderplaats and Yoshida 
(1985), and Fiacco and Ghaemi (1982). It is suspected that problems discussed by 
Robertson and Gabriele (1987) and Barthelemy and Sobieski (1983) also possess Case 3 
changes in the active set. Powell (1985) has studied the performance of RQP methods 
when the gradients of the constraints are linearly dependent. The test problems that were 
used by Powell can also be modified to be sensitivity test problems. 

We can expect to find other test problems (Case 3) when we begin to study more 
engineering test problems. Many engineering optimization problems are fully constrained 
at the optimum. When a new constraint becomes active for a fully constrained problem we 
will either have a Case 3 change in the active set or loose the feasible region(Case 4). 

To generate test problems for Case 3 changes in the active set the algorithm 
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described in section 6.4. 1 can be used with the following modification. In step 4 when the 
new constraint is added to the active set it will have to be linearly dependent with the 
constraints in the active set. To accomplish this the gradient of the new constraint at x(pl) 
can be constructed as a linear combination of the gradients of the other active constraints. It 
may be possible with further work to control which constraint in the active set is replaced 
when the active set changes. 

An alternative, less general procedure that can be used to create test problems with 
linearly dependent constraint gradients is illustrated for a two variable problem in figure 
6.7. Figure 6.7 (a) shows a simple two variable optimization problem. The problem has 
an elliptic objective function and is subject to an equality constraint that changes as the 
parameter p changes. There are also variable bounds present. For p = pO the optimum is 
located at the intersection of the equality constraint and the variable bound X2 ma x< When p 
is changed to p = pint the optimum is then located at the intersection of the equality 
constraint and both X2max and xi m i n . At this point the gradients of the constraints are 
linearly dependent which causes the linear independence assumption of the Kuhn-Tucker 
conditions not to hold. When p changes further to p = pi as shown in figure 6.7 (b) the 
optimum is now located at the intersection of the equality constraint and xi m i n . Thus as p 
moves from pO to pi there will be a change in the active set with the constraint gradients 
being linearly dependent. This procedure can be generalized to more dimensions. The 
discontinuity in df*/dp can be modified by varying the eccentricity of the ellipses. 



Figure 6.7 The Creation of Test Problems with Linear Dependencies 
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7 . 


Conclusions and Recommendations for Future Research 


7.1. CONCLUSIONS 

In this work we have proposed a new method for estimating parameter sensitivity 
based on the Recursive Quadratic Programming method. The new method approximates 
the sensitivities using a differencing formula and can be shown to be equivalent to a 
modified Kuhn-Tucker method. The method appears to be very competitive with existing 
methods when measured in terms of the number of function evaluations required to 
calculate a parameter sensitivity. It does not require the calculation of second order 
derivatives, but uses the BFS method or SRI method for developing an approximation to 
the Hessian of the Lagrangian. 

The choice of the variable metric update (BFS or SRI) effects the amount of work 
required to solve the perturbed problem, because different updates provide Hessian 
approximations of differing accuracies. Different variable metric updates also effect the 
speed with which the RQP algorithm can solve the problem. 

Initial testing of the algorithm against problems with known sensivities has shown 
that the method can adequately estimate the derivatives. The central difference 
approximation seems to provide the best results, particularly when the Hessian 
approximation is updated during the RQP iterations at the perturbed point. 

Using the RQP method to solve the optimization problem may be beneficial in many 
applications. This is because the RQP method has been shown to be one of the best 
general purpose algorithms for solving nonlinear programming problems. The use of 
variable metric updates allow the RQP method to solve problems where the Sequential 
Linear Programming (SLP) method fails. The RQP method can solve problems with very 
nonlinear constraints, and the RQP method performs the best when there are many active 
constraints at the optimum of the problem. 

Chapter 6 has discussed potential problems occurring when there are changes in the 

set. Chapter 6 also presented some techniques that can be used to deal with these cases. 

We observed discontinuities of the sensitivity derivatives in Chapters 5 and 6 when there 

were changes in the active set of constraints. If we are using sensitivity derivatives to make 

extrapolations, we now know Case 3 will cause the largest discontinuities (i.e. step 

discontinuities in the Lagrange multipliers) in the derivatives, and Case 4 will cause no 

solution to the proposed constraints to exist. We have also shown that the discontinuities 

in dx*/9p and 5u*/9p, that occur when the active set changes, make our prediction of 

which constraint will enter the active set second veiy difficult without reoptimizing the 
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problem. 

7.2. RECOMMENDATIONS FOR FUTURE RESEARCH 


The next step in the testing of the new RQP sensitivity algorithm will be to expand 
the test set to include a larger variety of test problems. We will need to include problems 
that have more variables, and also problems that demonstrate the behavior associated with 
Case 3 and Case 4 changes in the active set. We will also need to expand the test set to 
include engineering test problems for which the Hessian of the Lagrangian is not readily 
available. 

As we observed, the perturbation Ap can effect the accuracy of the derivatives that 
we calculated and we will have to investigate methods to improve the choice of Ap. In our 
initial testing we always let RQOPT use two iterations to solve the perturbed problem. 
Further tests are needed investigate the effectiveness of the algorithm when RQOPT is only 
allowed one iteration to solve the perturbed problem. 

We observed that the Hessian of the Lagrangian improved if we allowed the 
approximation to be updated during the reoptimization. More experiments need to be 
conducted to find how to best update the Hessian approximation. We also need to find 
when we should switch from using two iterations to solve the perturbed problem to using 
one iteration to solve the perturbed problem. Further work is also needed to identify when 
the Hessian approximation has converged. If the Hessian approximation converges to the 
true Hessian of the Lagrangian then the Hessian approximation may be used in the Kuhn- 
Tucker sensitivity equations, however 3V x L/dp will still need to be calculated. Since 
V X L * 0 the calculation of dV x L/9p may be subject to numerical noise and the RQP 
sensitivity algorithm may perform better than the Kuhn-Tucker method. 

The initial testing of the SRI update was very encouraging in terms of the 
convergence of the Hessian approximation to the True Hessian. However the SRI update 
proved to be unstable when used for test problem 4. We will need to investigate some 
method that can be used to stabilize the SRI update or find a set of rules that can be used 
that will only allow the SRI update to be used when the updated Hessian approximation 
will be stable. 

Using the RQP method to approximate the Hessian of the Lagrangian may be 

improved further if a Hybrid MOM/RQP algorithm is used to solve the original 

optimization problem. A Hybrid MOM/RQP algorithm would use the Method of 

Multipliers (or Augemented Lagrangian) algorithm for the first few stages, (build an 

approximation of the Hessian of the Lagrangian) then switch to using the RQP method. 

When the RQP method is used the approximation of the Hessian of the Lagrangian from 
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then MOM method can be used as the initial Hessian approximation, this should help the 
RQP method quickly solve the problem, and also obtain a good approximation of the 
Hessian of the Lagrangian. 

In summation we have seen that the new RQP sensitivity algorithm can find 
parameter sensitivities. We still need to investigate if this method will be superior to the 
Kuhn - Tucker method for general problems. Section 3.2.2 demonstrated that there will be 
a trade off with regarding the number of function evaluations required by the Kuhn-Tucker 
method versus the RQP method. We will have to conduct more experiments to study the 
general accuracy that can be expected from the RQP sensitivity algorithm. 

Using first order extrapolations can provide unsatisfactory results when the 
functions are nonlinear with respect to p. This situation is illustrated in Figure 7.1, which 
shows the effect of variations in p 3 on xj*. On this plot a linear and quadratic extrapolation 
are presented as well as the actual values of the optimum of xj*(p 3 ). It can be seen that the 
linear extrapolation does not provide a good estimate of the new value of xi* when p 3 
changes, however the quadratic approximation provides an accurate estimate of xi*(p3) up 
until the point where the active set changes. 

If good estimates of the second derivatives can be found then more accurate 
estimates of the behavior of the optimum can be made by using second order 
extrapolations. 

There are few available methods to calculate second derivatives of the optimum with 
respect to parameter variation but it is possible to predict d 2 f*/dpi 2 for some problems. 
However the only published algorithm that was found for calculating 9 2 x*/9pj 2 requires 
third order derivatives, which are seldom available for engineering problems. Thus, an 
algorithm based on the central differencing variants of the RQP sensitivity algorithm is 
proposed that can be used to calculate second derivatives. 


When using the central difference approximation an estimate of the second 
derivatives can be calculated from 


d 2 y(pj) y + - 2y + y ~ 
dpi 2 (Api)2 


(7.1) 


Where y can represent f*, x*, u*, or g*, and the estimates of f + >*, x + »‘, u + >", and g + >' 
calculated in steps 2 and 4 of the RQP sensitivity algorithm are substituted appropriately in 
(7.1). When using a central differencing approximation this option can be effective for 
indicating when curvature is present, but may not be able to accurately predict the true value 

of the second order information. It should be noted that this procedure may have the most 
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difficulty in predicting d 2 u(pi)/8pi 2 because the RQP algorithm produces more accurate 
estimates of f + >* and x + >- than u + >‘. 

Second derivatives might also be useful for predicting when constraints enter the 
active set. A prediction of when a constraint enters the active set can be identified using the 
following model of the behavior of the constraints that are not in the active set 

gnew = g(x*,p) + J^Api + j^-(Api) 2 (7.2) 



Figure 7.1 A Comparison of Quadratic to Linear Extrapolations for the Sensitivity 

Approimations. 


As a final note to this report we will discuss using the RQP method and the RQP 
sensitivity algorithm in multilevel decomposition algorithms. Because multilevel 
decomposition algorithms solve the same subproblems for different values of system level 
parameters, we expect that as the multilevel decomposition algorithm converges (after the 
proper active set has been identified) that the true Hessian of Lagrangian in the subsystem 
will not change very drastically. If the RQP method is used to solve the subsystem level 
optimizations and perform the sensitivity analysis, we expect that the Hessian 
approximation will converge to the true Hessian of the Lagrangian as the multilevel 
decomposition converges. If we use the approximation of the Hessian of the Lagrangian 
from the previous subsystem optimization as the initial Hessian approximation for the next 
iteration, a better approximation of the Hessian of the Lagrangian should be obtained when 
we solve the new subsystem problem. 

The above discussion implies that using the RQP method in conjunction with the 
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multilevel decomposition method may mean that on the first few system level iterations the 
sensitivities calculated at the subsystem level may be inaccurate, but the accuracy of the 
sensitivity derivatives will improve after each iteration at the system level. In our opinion 
this behavior is suitable for use with multilevel decomposition since far from the optimum it 
is often not advantageous (or necessary) to perform exact line searches or have exact values 
of the gradients. However, as the solution is approached, we need to obtain more accurate 
derivatives and perform more exact line searches to obtain acceptable convergence. 
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Appendix 1 

Test problems 


This section will present a discussion of the test problems that were used in the 
initial testing. Selected plots of the optimum sensitivity are provided to show how the 
optimum varies when the parameter is perturbed. The plots will also demonstrate how 
changes in the active set effect the optimum sensitivity. 

On each plot that is presented, the base value of the parameter will be indicated by a 
vertical line indicating the value at which the sensitivity analysis was performed. Active set 
changes will be indicated by a vertical line and a label indicating the constraint that enters or 
leaves the active set when the parameter is perturbed from the base point. On the plots of 
the optimum objective function vs the parameter, a linear extrapolation will be indicated by 
a line showing the predicted value of the objective function using equation 2.8. 

The plots of the behavior of the optimum can be used as a tool when diagnosing the 
behavior of various algorithms used to predict parameter sensitivities. Nonlinearity in the 
paths f*(x*,p), x*(p), and u*(p) can be seen in these plots, this nonlinearity can be used 
to explain why (or how far) linear extrapolations are valid. The discontinuities in the 
sensitivities that occur when the active set changes can also be seen and we can use these 
discontinuities to establish regions where the extrapolations are valid. 
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Problem 1 : 

Minimize f(x) = (xi - pi) 2 + (X2 - P3) 2 
subject to: gi(x) = 2xi - x^ - p2 > 0 

g2(x) = P3 * 0.8x j - 2x2 > 0 
Variable bounds [qJ < x < 

Starting Point for Optimization x° = (0,3) P® = (3,1,3) 
Optimum Point: f(x*(p 0 )) = 1.25 x*(p°) = ( 2.5, 2.0) 

Both constraints are active: u*(p°) = ( 0.3,0.4) 

Hessian of the Lagrangian H=£ 2 'q^ 2^] 


Sensitivity derivatives 



dx* _ 

3pT“ 

1 1 
O O 

b b 

8u* | 

^T = l 

varied pi from 1.5 to 4.8 

!r- 

8x* 

<5p2~ 

[-0.1 

Lo.2. 

1 9u * _ 
J ^p2 

r-0. 13041 

L 0.0008 J 

varied pi from -3.0 to 1.5 

It- 

dx* _ I 

1 

r 12 i 

Lo.6j 

9u* l 

cipi -- l 

' .4048 1 
.-0.5896J 

varied pi from 2.2 to 3.95 


Special features: Hessian matrix is available. Quadratic objective function and 

quadratic constraints. Active set changes are introduced for large changes in 
the parameters. Problem is fully constrained at the optimum. Plots of the 
behavior of the optimum are presented in figures A. 1 and A.2. 

Constraint gi leaves the active set when pi > 4.5, and constraint g2 leaves 
the active set when pi < 2.0. 

Constraint g2 leaves the active set when p3 > 3.888, and constraint gi 
leaves the active set when p3 < 2.4322. Since u*(p3) is nonlinear a linear 
prediction of when the constraint will leave the active set will not be accurate 
for large variations in the parameter. 
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Problem 2: 


J5 i i I 


- 5 ' 

Minimize f(x) = 0.5x7] 15 1 

x + x T 

5pi 

Ll 1 5 J 


L 5 J 

subject to: gi(x) = pixi + X 2 + X 3 - 

P2 — 0 



g 2 (x) = xi + 2 x 2 + 3 x 3 - 4.7 - pi > 0 



r-iOi 


r 20 i 

Variable bounds 

-10 

L-10J 

< x < 

20 

L20J 


Starting Point for Optimization x° = (1.1, 1.2, 1.3) p° = (1,3) 
Optimum Point: f(x*(p°)) = 25.5 x*(p°) = ( 1.0, 1.0, 1.0) 
gl active at the optimum: u*(p°) = ( 12 . 0 , 0 . 0 ) 


Hessian of the Lagrangian H 



1 

5 

1 


1 I 
1 

5 J 


Sensitivity derivatives 


df_ 

dpi 


-7.0 


9x* 

3pT 


r 2.0833333 -i 
-2.1666666 
L -.9166666 J 


3u* _ r -4.666666661 

3pT“- L 0.0 J 


varied pi from 0.8 to 1.2 


df 

dp2 


12.0 


9x* 

dp2 


0.333333 I 
0.333333 
L0.3333333J 


9u* _ r2.333333~| 
<?P2 L 0 J 


varied p 2 from 2.7 to 3.2 


Special features: Hessian matrix is available. Quadratic objective function and 
linear constraints. Parameter 2 is a right hand side perturbation. 

A plot of the behavior of the optimum with respect to pi is presented in 
figure A.3. For pi > 1.04889, constraint g 2 enters the active set. We can 
see a change in sign of d^f*/dpj from a region of positive curvature to a 

region of negative curvature. We can also see the discontinuity of the slope 
of 9x*/9p and du*/ 9 p when the active set changes. 

A plot of the behavior of the optimum with respect to p 2 is presented in 
figure A.4. For pi < 0.95, constraint g 2 enters the active set.. Again we 
can see the discontinuity in 9 x*/ 9 p 2 and 9 u*/ 9 p 2 . 


C-3- 
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Problem 3 (Common name Rosen and Suzuki Test Problem): 
Minimize f(x) = x^ + x^ + 2 x^ + x^ - 5 xi - 5x2 - 21x3 + 7x4 + 100 

2 2 2 2 

subject to: gi(x) = (-Xj - x 2 - x 3 - x 4 - xi + X2 - X3 + X4 )/8 + pi > 0 
g 2 (x) = (- x i * 2 x 2 - x^ - 2 x^ + xi + X 4)/10 + P2 ^ 0 
g3( x ) = (-2xj - x^ - X3 - 2xi + X2 + X4)/5 + p3 > 0 


Variable bounds 


r-10i 


[20-1 

-10 

<x < 

20 

-10 

20 

L-10J 


L20J 


Starting Point for Optimization x° = (0.0,0.0,0.0,0.0) p° = (1,1,1) 
Optimum Point: f(x*(p 0 )) = 56.0 x*(p°) = ( 0.0, 1.0, 2.0, -1.0) 


gl and g3 active: 

Hessian of the Lagrangian H=) 


u*(p0) = (8.0, 0.0, 10.0) 


[-12 0 0 0 
0 8 0 0 
0 0 10 0 
L0 0 0 4 


df 

dpi 


df 

dP2 


df 

dpi 


Sensitivity derivatives 
dx* 


= - 8.0 


= 0.0 


3pT : 


dx* 

<5p2" 


r -1.14719841 
0.33428030 
-0.2279022 
L -3.5403608 J 


ro.o 

0.0 

0.0 

L0.0J 


= - 10.0 


dx* 

<5pT 


du* 

<5pT" 


du* 

<5i>2' 


1.3057920] 

0.5460588 

du* 

-1.0541310 

dpi ~ 

L-2.3741690J 


r-67.3428300"| 

0.0 

55.4605800 


0.0i 

0.0 

L0.0. 


-67.3428300- 

0.0 

. 55.4605800 J 


varied pi from 0.86 to 1.12 


varied p 2 from 0.85 to 1.15 


varied pi from 0.86 to 1.12 


Special features: Hessian matrix is available. Quadratic objective function and 

quadratic constraints. Problem has been used in the other studies of optimal 
design, the parameters that we are studying are right hand side 
perturbations. 


94 



Figure A.5 and A.6 present plots of the variation of the optimum with 
respect to pj. The following changes in the active set are introduced: when 
pi < .855, constraint g3 leaves the active set; when pi > 1.057, constraint 
g2 enters the active set; when pi > 1.078, constraint gi leaves the active set. 
The plot shows that f*(pi) is nonlinear. It is difficult to see the 
discontinuity in dx*/9pi but we can clearly see the discontinuity in 9u*/9pi. 
Figure A.7 presents a close up view of the behavior of X3 and X4. We can 
see that X3 is a nonlinear and a piecewise continuous function of pi. The 
discontinuities take place when the active set changes. With the resolution 
in figure A.5 this behavior was very difficult to see, however with the 
enlarged view this becomes easy to see. 

A plot of the behavior of the optimum with respect to p2 is presented in 
figure A. 8. When p2 < 0.9, constraint g2 enters the active set. In figure 
A. 8 we can see a sharp discontinuity in du/dp and we also can see that large 
errors would be introduced if the value of the objective function was 
extrapolated from the base point to after the active set changed. 



Problem 4: 


5 5 5 5 3 

Minimize f(x) = Xe;xj + Y £cii x i x i+ £djx? 

j=l i=lj=l i=l Jj 


5 

subject to: gi(x) = £ayxi - bi > 0 i=l,10 

j=l 


Where the values of ay, bi, cy, di, ej are constants that can be found in (Coville 
1969, Himmelblau 1972, Eason and Fenton 1974, or Sandgren 1977) 


Variable bounds 


r°i 


r20-| 

0 


20 

0 

< X < 

20 

0 


20 

LoJ 


l-20-l 


Starting Point for Optimization x° = (0.,0.,0.,0.,1.) p° = (b 3 , bs, b6, b9)=(-.25,-4.,-l.,5.) 
Optimum Point: f(x*(p°)) = -3.234866708 

x*(p0) = (0.3, 0.3334676, 0.4,0.4283101,0.2239649) 


g3> g5> g6> and g 9 are active: 

u*(p°) = (0.,0.,0.5 17404, 0.,0.306111,1. 183954594, 0.,0., 0.010390,0.) 


Hessian of the Lagrangian 



r6.72 

-4.0 

-2.0 

6.4 

-2.0 


-4.0 < 

?.4006 

-1.2 

-6.2 

6.4 

H= 

-2.0 

-1.2 

4.4 

-1.2 

-2.0 


6.4 

-6.2 

-1.2 9.3418 

-4.0 


L-2.0 

6.4 

-2.0 

-4.0 

6.2688 


Sensitivity derivatives 


f -0.4 n 


df 

dpi 


= .517404 



0.097300 

- 0.2 

0.286025 
L- 0.067740-1 


r 0.407099 -i 
0u*_ -0.056624 
0.347303 
L 0.019084 J 


varied pi from -0.39 to -.019 


96 



df 

dp2 


= .306111 


dx* 

<5p2 


0.0 "I 
-0.147119 
- 0.0 

-.049166 
0.09817 J 


du* 

3i>2 


-0.056624-1 

.0505155 

-.061564 

L0.0024016J 


varied p 2 from -4.5 to -3.5 


df 

dp3 


1.18395 


3x* 

3p7 


-0.2 “I 
0.094082 
-0.35 
0.228817 
1-0.02931 34 


r .34303 *i 
du*_ -0.061564 
3p7“ 0.720695 
L 0.1 59647 J 


varied P 3 from - 1.2 to - 0.8 


^=.010390 

dp4 


dx* 

^P4 


0.0 -I 
-0.039759 
0.0 

0.0761408 
L. 15499 104-1 


rO.O 19085-1 
du*_ 0.002402 
Tfc- 0.159647 
L0.083860J 


varied p 4 from 4.0 to 6.0 


Special features: Hessian matrix is available. Cubic objective function and linear 
constraints. Active set changes are introduced for changes in P 4 . The 
parameters being studied are right hand side perturbation. This problem has also 
been studied by Fiacco et. al. (1974,1983). 

A plot of the sensitivity of the optimum with respect to pi is presented in figure 
A.9. Plots of the optimum sensitivity with respect to p 2 and p 3 are similar to 
those to pi. For variation in the first three parameters, the optimum objective 
function behaves linearly and there are no changes in the active set for the range 
of parameter perturbation that was studied. 

A plot of the sensitivity of the optimum with respect to p 4 is presented in figures 
A.lOandA.ll. Constraint g 9 leaves the active set when p 4 < 4.876. After 
constraint g 9 has left the active set further variation of p 4 has no effect on the 
optimum. For increasing values of p 4 ) errors would be introduced if a linear 
extrapolation was used. 
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«x*) 


Figure A. 5 Sensitivity with respect to p(l) for Problem 3 
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Lagrange Multipliers 


Figure A.6 Sensitivity with respect to p(l) for Problem 3 (cont) 
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Figure A.7 Plots of the Sensitivity of x(3) and x(4) for problem 3 p(l) 
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p(2) with p(l)=p(3)=l p(2) with p( I )=p(3)= 1 



safqmjBA uSissq 


p(2) with p(1)=p(3)=l P< 2 ) wilh P(')=P(3)=I 









Problem 4: 


5 5 5 5 - 

Minimize f(x) = £ejxj + £ £cjjXjXj+ £djx? 

j=l i=l j=l i=l 


5 

subject to: gi(x)= £aijXi-bi>0 i=l,10 
j=l 


Where the values of ajj, bj, cjj, dj, ej are constants that can be found in (Coville 
1969, Himmelblau 1972, Eason and Fenton 1974, or Sandgren 1977) 


Variable bounds 


roi 


r20-| 

0 


20 

0 

< X < 

20 

0 


20 

LoJ 


L20J 


Starting Point for Optimization x° = (0.,0.,0.,0.,1.) p° = (b 3 , bs, b6, b9)=(-.25,-4.,-l.,5.) 
Optimum Point: f(x*(p°)) = -3.234866708 

x*(p0) = (0.3, 0.3334676, 0.4,0.4283101,0.2239649) 
g3> g5» g6» and g 9 are active: 

u*(p0) = (0.,0., 0.517404,0., 0.3061 11,1. 183954594, 0.,0.,0.0 10390,0.) 


Hessian of the Lagrangian 


f6.72 -4.0 -2.0 6.4 
-4.0 9.4006 -1.2 -6.2 
H= -2.0 -1.2 4.4 -1.2 

6.4 -6.2 -1.2 9.3418 

L-2.0 6.4 -2.0 -4.0 


- 2.0 1 
6.4 
- 2.0 
-4.0 
6.2688-1 


Sensitivity derivatives 


df 

dpi 


= .517404 



-0.4 “I 
0.097300 
- 0.2 

0.286025 

L-0.067740J 


r 0.407099 1 
3u* -0.056624 

X7= 0.347303 
L 0.019084 J 


varied pi from -0.39 to -.019 
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df 

dpi 


= .306111 


ax* 

3i>2 


0.0 -I 

-0.147119 

- 0.0 

-.049166 
0.09817 J 


9u* 

^P2 


-0.056624-1 

.0505155 

-.061564 

L0.0024016J 


varied p 2 from -4.5 to -3.5 


df 

dp3 


1.18395 


r -0.2 l 


0X* 

^pT 


0.094082 

-0.35 

0.228817 

L0.029313J 


r .34303 -| 
3u*_ -0.061564 
3p7~ 0.720695 
L 0.159647 J 


varied P 3 from - 1.2 to - 0.8 


df 

dp4 


= .010390 


3x* 

ap4 


o.o -| 

-0.039759 

0.0 

0.0761408 
L. 154991 04 J 


r0.019085-i 
au*_ 0.002402 
0.159647 
L0.083860J 


varied p 4 from 4.0 to 6.0 


Special features: Hessian matrix is available. Cubic objective function and linear 
constraints. Active set changes are introduced for changes in P 4 . The 
parameters being studied are right hand side perturbation. This problem has also 
been studied by Fiacco et. al. (1974,1983). 

A plot of the sensitivity of the optimum with respect to pi is presented in figure 
A.9. Plots of the optimum sensitivity with respect to p 2 and p 3 are similar to 
those to pi. For variation in the first three parameters, the optimum objective 
function behaves linearly and there are no changes in the active set for the range 
of parameter perturbation that was studied. 

A plot of the sensitivity of the optimum with respect to p 4 is presented in figures 
A.lOandA.ll. Constraint g 9 leaves the active set when p 4 < 4.876. After 
constraint g 9 has left the active set further variation of p 4 has no effect on the 
optimum. For increasing values of p 4 > errors would be introduced if a linear 
extrapolation was used. 
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1. a grange Multipliers Design Variables l(x*) 












Ucsiun Variables 


Figure A. 10 Sensitivity with respect to P(4) for Problem 4 



p(4)=b(9) with p(l)=-.25, p(2)=-4, p(3)=-l 



p(4)=b(9) with p(l)=-.25, p(2)=4, p(3)=-l 
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Figure A. 11 Sensitivity with respect to P(4) for Problem 4 (corn) 
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Appendix 2 


THE ROCRE AND ROSEN SYSTEMS 

This section will discuss the system that was created for studying parameter 
sensitivities. The first subsection introduces the support program RQCRE that was written 
to simplify the construction of test problems. The second section will discuss the RQSEN 
system. The RQSEN system is a interactive program that acts as a post 
processor/sensitivity analysis module for the RQOPT program. 

RQCRE and RQSEN were set up to have some user friendly features. The 
programs have some user-freindly features however the programs will crash if, the user 
enters numbers in response to a promt for an alpha data type, and may crash if the user 
enters letters when the system is requesting numbers. 

A2.1 The RQCRE Support System 

The RQCRE program was written to reduce the time required to implement test 
problems. The RQSEN program requires approximately 30 arrays and a complicated main 
program to be written by the user. The RQCRE program automatically dimensions the 
proper arrays and writes the required calling programs. Using the RQCRE program 
reduced the time required to implement test programs during our initial testing. 

The RQCRE program is essentially a program that writes another program. The 
main features of the RQCRE program are; 

1. The program can be used in an interactive mode. 

2. The program writes the main calling program. 

3. The program can write an outline of the function subprogram. 

4. The program can be used to update the problem formulation. 

A structure chart of the RQCRE program (in CMS) is presented in Figure A2.1. 

The basic functions of the program modules are; 

RQCRE.EXEC - this module connects the proper output files to the proper unit 
numbers. 

RQCRE - is a FORTRAN program that can be used interactively to create a problem 
for submission to the RQSEN system. The input to RQCRE can either 
come from a data file or from the user. The output from RQCRE is a data 
file "data" that can be used as an input file for RQSEN, a MAIN FORTRAN 
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program ready for compilation, and a shell for the function subprogram. 
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Figure A2.1 A Structure Chart for the RQCRE Program. 

3 =N 

0 =LFSW 

0 =NEQ 

0 =LEQ 

2 =NINEQ 

2 =LINEQ 

2 =NPARM 

15 =MAXIT 

1 =1 START 

1 =IFIN 
50 =LPP 

6 =NOUTFL 

5 =ITR 

l.E-4 =PMIN 
l.E-8 =CAPFMN 

0 . 5 =DELTA 

10.0 =R 
.0001 =GAMMA 
1.0E-8 =EPSQP 
5.E-4 =EPSGRD 

2 =XDIF 

0 =NSCALE 

1. E-11 =ZEROM 

1.1 =X(1) 

1.2 =X (2) 

1.3 =X (3) 

-10. =XMIN (1) 

-10. =XMIN (2) 

-10. =XMIN(3) 

20.0 =XMAX (1) 

20.0 =XMAX(2) 

20.0 =XMAX (3) 

1.0 =PAFM(1) 

3.0 =PAFM(2) 

Figure A2.2 A Sample of a Data file (Test Problem 2) for the RQSEN Program. 
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data - the input/output data file that contains the algorithm parameters, starting point 
and initial values of the parameters for RQSEN. 

MAIN - FORTRAN program used as the main calling program when running 
RQSEN. 

FSUBI - a FORTRAN function that the user is required to modify by adding the 
definitions of the constraints and objective function. 

A sample of the data file for test problem 2 is presented in Figure A2.2. This file 
was written by the RQCRE program. This data file is used as an input to the RQSEN 
system to provide the programs with the values of the algorithm parameters, design 
variables, and initial values of the parameters. 

A sample of the a program written by the RQCRE program is provided in Figure 
A2.3. The program represents an implementation for test problem 2 (described in appendix 
1). The only modifications that were made to the program are, the objective function and 
the constraint definitions that were added to the code generated by RQCRE. 

A2.2 The RQSEN program 

This section describes the implementation of the RQSEN system. The first topic to 
be discussed is the capabilities of the RQSEN system. Next a description of the 
implementation is provided. The final topic presented in this section is a sample session 
from the RQSEN system. 

The basic capabilities of the RQSEN system are; 

1. The program can be used to solve optimal design problems. 

2. The program can be used to conduct convergence studies for various versions of 
RQOPT. 

3. The program can be used to calculate parameter sensitivity derivatives. 

4. The program can be used to conduct studies of large variations in the 
parameters. 

5. The program can be used to create parameter sensitivity plots that can be used, 
for trade off studies. 

The RQSEN system is currently implemented on the following systems, an IBM 
4341 under the CMS operating system and a microVAX under the VMS operating system. 
The programs are written in FORTRAN 77 and implemented in double precision. 

Figure A2.4 presents a structure chart for the RQSEN program (CMS 
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PROGRAM RTS02 


C 

C TEST PROBLEM 2 PH. D., CREATED BY TODD J. BELTRACCHI TO EXAMINE 
C CONVERGENCE OF THE HESSIAN APPROX AND CHANGES IN THE ACTIVE SET 
C THIS HAS A QUADRATIC OBJECTIVE FUNCTION AND LINEAR CONSTRAINTS 
C 

IMPLICIT REAL* 8 (A-H f O-Z) 

DIMENSION X (3) ,XMIN(3) ,XMAX (3) , SCALE (3) , H (3, 3) , DELFHG (3, 3) , 

1 FUNCT (3) ,P(3) ,U(8) ,V(1) ,XS (3) , DFHGS (3, 3) , FUNCTS (3) , PS (3) , US (8) , 

2 VS (1) ,DFDP(2) ,FUNCTP(3,2) ,DXDP(3,2) ,DUDP(8,2) ,DVDP(1,2) ,DFDPE(2) , 

3 DGDP (2,2) 

LOGICAL YESNO, ACTPT ( 8 ) , ACTPTS ( 8 ) 

CHARACTER* 3 YSN 
CHARACTER* 7 FILENM 
COMMON /OPTDAT/ D(400) 

COMMON /BFSWK/ DD(20) 

COMMON /PMINI/ PMINI (3) 

COMMON /PMAXI/ PMAXI (3) 

COMMON /PARMS/ PARM(2) 

INCLUDE (RQS) 

END 

Q ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 

FUNCTION FSUBI (X, IEVAL) 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION X (3) 

COMMON /NFEVAL/ NCE,NFE 
COMMON /PARMS/ P(2) 

IF (IEVAL .GT. 1) GOTO 2 

1 NFE=NFE+1 

C PLACE OBJECTIVE FUNCTION DEFINITION HERE 

FSUBI=2 . 5* (X (1) **2+X (2) **2+X (3) **2) +X (1) *X (2) +X (1) *X (3) +X (2) *X (3) 

1 +5. * (X (1) +P (1) *X (2) +X (3) ) 

RETURN 

2 NCE=NCE+1 

C PLACE LINEAR EQUALITY CONSTRAINTS HERE 
GOTO (10, 20, ) IEVAL 

10 FSUBI=P (1) *X (1) +X (2) +X (3) -P (2) 

RETURN 

20 FSUBI=X (1) +2 . *X (2) +3 . *X (3) -4 . 7-P (1) 

RETURN 

END 

Figure A2.3 A Sample Program (Test Problem 2) for RQSEN. 
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implementation). A brief explanation of each of the program modules is provided 

RQSEN.EXEC - this module connects the proper files to the proper unit numbers 
and prompts the user for the name of the program to be run. 

MAIN - the main calling program. This module calls RQOPT and RQSEN, this 
module also reads in the starting point and algorithm parameters for RQOPT 
and allows the user to save solution point to a data file for later use, i.e. a 
sensitivity study at a later time. 

FSUBI - the function subprogram that defines the objective function and 
constraints. 

data - the data file that contains the algorithm parameters and values of the design 
variables. 

RQOPT - implementation of the RQP method described in Beltracchi and Gabriele 
(1987 b). 

RQSEN - the main driving routine for a sensitivity analysis. 

WRTPT - A utility program for writing the design point, Lagrange multipliers and 
values of the objective function and constraints to a summary file. The 
summary file is in a form which can be read by a plotting program to 
graphically display the sensitivity information. 

solf - a data file used to communicate optimum design points to a plotting program. 

PARTP - a subroutine that uses either forward, central or user supplied routines to 
calculate the partials with respect to the design parameters of the objective 
function and constraints. 

PARSEN - a subroutine that calculates 9x*/9p, 9u/9p, df/dp and dg/dp by either 
forward or central differencing. The user can specify the perturbation size 
and the number of iterations that RQOPT uses to solve the perturbed 
problem. 

PARSEN2 - a subroutine similar to PARSEN but this imnplements modified central 
differencing. 

RQDR - is the routine that controls the execution of RQOPT for reoptimization. 

PARVAR - the subroutine used to conduct studies of large variations in problem 

parameters, if parameter sensitivity derivatives are available then the location 
of the starting point for the reoptimization is approximated by. 

^initial = *(P) + ^ Ap 

PRVRDR - is the routine that controls the execution of RQOPT for reoptimization. 
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The rest of this section describes the steps involved in using the RQSEN system. 

The first step is for the user to create the necessary FORTRAN code to define the 
main calling program and the function subprogram (similar to the one in Figure A2.3). The 
second step is to set up a data file (similar to the one presented in Figure A2.2) with the 
values of the design variables and algorithm parameters. These first two steps can be 
performed with the aid of the RQCRE preprocessor. 

Once the calling program, function subprogram, and the data file are defined, the 
user can run the RQSEN system to conduct a study of the sensitivity of the optimum of the 
problem or to study the convergence of the problem. 

A sample of some convergence plots are presented in Figure A2.5 & 6. These plots 
can be used to assess the convergence criteria of various algorithms, for example Figure 
A2.5 shows that the BFS version was able to solve the test problem faster than the version 
that used the SRI update. However the SRI update found a more accurate estimate of the 
optimum, and once the region of the minimum was located the convergence for the SRI 
update was better than the BFS update. 

The best way to explain how the program can be used to conduct parameter 
sensitivity studies is to provide a sample session (see Figure A2.7) that was run under the 
IBM version of RQSEN. The next several paragraphs describe the output in Figure A2.7, 
all user responses are shown in italics. 

The sample problem 2, described in appendix 1, is solved. The modified SRI 
update is used to approximate the Hessian of the Lagrangian, and modified central 
differencing is used to calculate the sensitivity derivatives. The first step for running 
RQSEN to invoke the exec file to start the program, this is done by entering rqsri 4. The 
first prompt asks the user for a name of a data file to store the results in the user responds 
with rts02. This data file can be used to maintain a summary of all optimum points. The 
next prompt is for the name of the program to be run, in our example program rts02, an 
implementation of test problem 2 (see appendix 1) was run. The program and data file 
were presented in Figure A2.3 and A2.2. 

When the program begins the first prompt is for the name of the file with the input 
data. A response of rts02 is entered, the program then reads in the data from file r ts02. 
The next prompt asks if there is an approximation to several arrays available. These arrays 
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Iteration 

Figure A2.5 A Sample Plot comparing the Convergence Characteristics of the 

BFS and SRI Updates. 



Iteration 

Figure A2.6 A Sample Plot Showing the Convergence of the Design Variables. 
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are saved if the optimum has already been found, a response of "n " for no is entered 1 . 
Now the program invokes the RQOPT program to locate the optimum of the problem. A 
summary of the output from RQOPT (during the optimization) is presented in Figure 
A2.7(a-c). After the problem has been solved some final statistics from RQOPT are 
presented along with an approximation to the Hessian of the Lagrangian, in both the LDL T 
format and the unfactored form. These are shown in Figure A2.7(b & c). 

After the problem has been solved by RQOPT control is returned to the 
preprocessor (see Figure A2.7(c)). The preprocessor provides the user with a choice of 
being able to save the optimal point. In the example the response was "y " for yes was 
entered, next the user is asked if he wants to save the final point in the same file as the 
initial point, a response of "n " for no was entered 2 . Next the user is asked to supply a new 
name of the data file to store the point in, a response of rts02s was entered. The data file 
was then written and the user asked if they wanted the gradients and Hessian 
approximation to be written to the file, a response of "y " was entered 3 . The preprocessor 
next asks if the user wants to perform a sensitivity analysis, a response of "y " for yes is 
entered. 

Now control is passed the the RQSEN program (see Figure A2.7(c-e)). The first 
question asked by RQSEN is if the user wants the solution points written to a solution file a 
response of "n " is entered. The next question asked is for a value of epsilon to calculate 
the partial derivatives of the problem functions. RQSEN then calculates the partials of the 
objective function and constraints, the derivative of the objective function is then calculated 
by equation 1.20. 

The next step in the sample output is the calculation of the partials of the optimum 
design variables and Lagrange multipliers with respect to the first parameter for the 
problem. Again the user can specify the size of the perturbation of the parameter and the 
number of iteration that RQOPT is allowed to use for solving the perturbed problem. In 
this example central differencing is used, and Hessian updating is allowed, notice that in 
Figure A2.7(d) the Hessian approximation has converged. 

Figure A2.7(e) shows the values of the parameter sensitivity derivatives that were 
calculated by RQSEN. The gradient of the objective function df*/dp was calculated by 3 

iRQCRE or RQSEN will accept either "YES", "YE", "Y", "yes", "ye", or "y" for a yes responce and 
"NO", "N", "no", "n”, null for a no responce. 

2 If the user responded yes then the final point would overwrite the initial point in the data file. 

3 The Data file can now be used as an input to RQSEN for a sensitivity analysis performed using the 
gradients and Hessian approximation that were found when solving the problem. 
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different approximtions, all values are reported for comparison. The next option of the 
program is to study finite perturbations in the parameter. This option can be used to 
calculate optimal designs for different values of the parameters and then to write the optimal 
design points to a data file (see Figure A2.7(e&f)). The rest of the sample output involves 
terminating the program. 

When the user saves the optimal points plots of the optimum sensitivity can be 
made. An example of the sensitivity for problem 2 when pi is varied is presented in Figure 
A2.8. 




Figure 3.8 A Plot of the Sensitivity of the Optimum Objective Function and Optimum value 
of the Design Variables for Problem 2 when pi is perturbed 


117 




rqsrl 4 

INPUT NAME OF THE FILE TO HOLD THE SENSITIVITY DATA 
r ts02 

FILEDEF 9 DISK RTS 02 DATA A1 ( LRECL 800 

INPUT NAME OF THE PROGRAM TO BE RUN BY RQSEN 
rts02 

LOAD RTS 02 SRI 4 RQOPT OPTQP RQSEN PARSEN2 PARVAR WRTPT ( CLEAR START 
Execution begins . . . 

★ ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 

* ENTERING RQOPT/RQSEN PREPROCESSOR * 

★ ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★•A-* 


INPUT NAME OF THE FILE WITH INPUT DATA (XXXXXXX) 
rts02 

OPENING FILE, RTS 02 ON UNIT= 3 FOR INPUT OF PROBLEM DATA 

ARE FUNCT, H, U, V,ACTPT, DELFHG IN FILE RTS02 (YES/NO) 
n 

RQOPT VERSION 2.12 - EXPERMENTAL 2/9/88 
STARTING INFORMATION 


INPUT PARAMETERS: 

NUMBER OF VARAIBLES = 3 

OPTIMIZATION TO BE PERFORMED ON A NON-LINEAR OBJECTIVE FUNCTION 
NUMBER OF EQUALITY CONSTRAINTS = 0 

NUMBER OF LINEAR EQ CONSTRAINTS = 0 

NUMBER OF INEQUALITY CONSTRAINTS = 2 

NUMBER OF LINEAR INEQ CONSTRAINTS = 2 

MAXIMUM NUMBER OF ITERATIONS =100 
LINES PER PAGE = 50 

ITR = 1 

NUMBER OF OUTPUT FILE = 6 


DELTA 

R 

GAMMA 

EPSILON FOR THE QP 
EPSILON FOR THE GRAD 
DIFFERNCING TYPE 
INITIAL VAULE OF CAPF 
MINIMUM VALUE OF CAPF 
MINIMUM NORM OF P VECTOR 
SCALING PARAMETER NSCALE 
NUMBER OF PARAMETERS NPARM 
INDEX PARAMETER VALUE 

1 1.00000000000000000 

2 3.00000000000000000 


0 . 500000E+00 
0 . 100000E+02 
0.100000E-04 
0.100000E-09 
0 . 100000E-03 
2 

0.000000E+00 
0 . 100000E-06 
0 . 100000E-03 
0 
2 


I XMIN 

1 0 . 000000E+00 

2 0 . 000000E+00 

3 0 . 000000E+00 


XMAX 

0.100000E+03 
0 . 100000E+03 
0 . 100000E+03 


SCALE 

0 . 100000E+01. 
0 . 100000E+01 
0 . 100000E+01 


Figure A2.7(a) A Sample of RQSEN for Test Problem 2. 


118 



THE HESSIAN APPROXIMATION IN LDL (T) FORM 
ROW 1 0.100000E+01 

ROW 2 0. 000000E+00 0.100000E+01 

ROW 3 0.000000E+00 0.000000E+00 0.100000E+01 


ENTERING DRIVER ROUTINE 

ITERATION 0 : 0 

OBJECTIVE FUNCTION = 0 . 33160000E+02 FUNCTION EVALUATIONS= 1 

CONSTRAINT EVALUATIONS= 14 

H (I) = ?0 V(I) 


PAGE 1 


G (I) >= ?0 
0.600000E+00 
0 . 170000E+01 


U(I) 

0 . 000000E+00 
0 . 000000E+00 


INDEX X(I) 

1 0 . 1100000E+01 

2 0 . 1200000E+01 

3 0 . 1300000E+01 

CHECKING CONVERGENCE THE NORM OF P= 1.96446237582566496 

ITERATION 1 : 0 

PAGE 1 OBJECTIVE FUNCTION = 0 . 25535366E+02 FUNCTION EVALUATIONS= 8 

CONSTRAINT EVALUATIONS= 19 

H (I) = ?0 V(I) 


G (I) >= ?0 
A-.288658E-14 
A0.488060E+00 


U(I) 

0 . 000000E+00 
0 . 000000E+00 


INDEX X(I) 

1 0 . 9059701E+00 

2 0 . 1000000E+01 

3 0 . 1094030E+01 

CHECKING CONVERGENCE THE NORM OF P= 0.345110459418981358 

ITERATION 2: 1 

PAGE 1 OBJECTIVE FUNCTION = 0 . 25503133E+02 FUNCTION EVALUATIONS= 16 

CONSTRAINT EVALUATIONS= 23 

H (I ) = ?0 V(I) 


G (I) >= ?0 
A0.577316E-14 
A0.244030E+00 


U(I) 

0 . 586085E+01 
0 . 659360E-01 


INDEX X(I) 

1 0 . 1027985E+01 

2 0 . 1000000E+01 

3 0 . 9720149E+00 

CHECKING CONVERGENCE THE NORM OF P= 0 . 395769395140752175E-01 

ITERATION 3 : 2 

PAGE 1 OBJECTIVE FUNCTION = 0 . 25500000E+02 FUNCTION EVALUATIONS= 23 

CONSTRAINT EVALUATIONS= 25 

H (I ) = ?0 V(I) 


G (I) >= ?0 U(I) 
A0.288658E-14 0.120000E+02 
A0.300000E+00 -. 288658E-14 


INDEX X(I) 

1 0 . 1000000E+01 

2 0 . 1000000E+01 

3 0 . 1000000E+01 

CHECKING CONVERGENCE THE NORM OF P= 0 . 395476598531791093E-08 
CONVERGENCE ACHIEVED 


FINAL STATISTICS 
CONVERGENCE ACHEIVED 

IN 29 FUNCTION EVALUATIONS 

4 FUNCTION GRADIENTS 
25 CONSTRAINT EVALUATIONS 

2 CONSTRAINT GRADIENTS 

3 I ITERATIONS 

WITH 0.00000000E+00 BEING THE MAXIMUM CONSTRAINT VIOLATION 

Figure A2.7(b) A sample of RQSEN for test problem 2. 
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THE HESSIAN APPROXIMATION IN LDL (T) FORM 
ROW 1 0 . 450000E+01 

ROW 2 0 . 444444E+00 0.211111E+01 

ROW 3 0 . 111111E+00 0 . 842105E+00 0.294737E+01 

THE HESSIAN APPROXIMATION UPPER TRIANGLE 
ROW 1 0 . 450000E+01 0.200000E+01 0.500000E+00 

ROW 2 0 . 300000E+01 0.200000E+01 

ROW 3 0 . 450000E+01 

PAGE 1 OBJECTIVE FUNCTION = 0 . 25500000E+02 FUNCTION EVALUATIONS= 29 

CONSTRAINT EVALUATIONS= 25 
INDEX X(I) H (I) = ?0 V(I) G (I) >= ?0 U(I) 

1 0 . 1000000E+01 A0.288658E-14 0.120000E+02 

2 0 . 1000000E+01 A0.300000E+00 O.OOOOOOE+OO 

3 0 . 1000000E+01 

DO YOU WISH TO SAVE THE FINAL DATA (YES/NO) ? 

y 

DO YOU WISH TO USE THE SAME FILE RTS02 (YES/NO)? 
n 

INPUT NAME OF THE FILE FOR STORAGE OF DATA (XXXXXXX) 

rts02s 

OPENING FILE, RTS02S ON UNIT= 4 TO STORE PROBLEM DATA 

DO YOU WANT FUNCT, H, U, V, DELFHG, ACTPT WRITTEN TO RTS02S (YES/NO)? 

y 

DO YOU WANT TO PERFORM A SENSITIVITY ANALYSIS (YES/NO) ? 

y 

WELCOME TO RQSEN1 . 0 

A SENSITIVITY ANALYSIS PROGRAM FOR RQOPT 
LAST MODIFIED APRIL 28 1988 

DO YOU WANT TO WRITE THE SOLUTION POINTS TO FILE= 9 (YES /NO) ? 

n 


THE DERIVATIVE OF THE OBJECTIVE FUNCTION WITH 
RESPECT TO ALL PARAMETERS WILL BE CALCULATED 

INPUT EPSP FOR THE CALCULATION OF DF/DP 
2 

0001 


THE DERIVATIVE OF OBJECTIVE FUNCTION W.R.T. P 
-0 . 6999999989E+01 0 . 1200000000E+02 

DO YOU WANT TO STUDY FINITE PERTURBATIONS 

ENTER PARAMETER NUMBER OR (-1 OR CTRL Z) TO CALCULATE GRADS? 
2 

-1 

DO YOU WISH TO FIND PARTIALS OF THE DESIGN VARIABLES AND 
LAGRANGE MULTIPLIERS W.R.T. PARAMETER (NUMBER OR -1 TO END)? 

2 

1 


Figure A2.7(c) A Sample of RQSEN for Test Problem 2. 
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ENTER EPSP FOR THE GRADIENT CALCULATION? 

? 

.0001 

PERFORMING A SENSITIVITY ANALYSIS FOR PARM ( 1) 

ASSUMING BASE POINT IS STABLE 

BASE POINT VALUE PARM ( 1 ) =0 . 1000000000000000E+01 

PERTURBED VALUE OF PARM ( 1) =0 . 1000100000000000E+01 

ENTER THE NUMBER OF ITERATIONS FOR RQOPT? 

? 

2 

★ ★★★★★★★★★★★★★★★★★★★★★★★★★ ENTERING RQOPT ^^^^^t******'**'*''*''*''*'**'* - '*''* 


ENTERING DRIVER ROUTINE 

ITERATION 0: 0 PARAMETER ( 1)= 0 . 10001000E+01 

PAGE 1 OBJECTIVE FUNCTION = 0 . 25500500E+02 FUNCTION EVALUATIONS= 7 

CONSTRAINT EVALUATIONS= 26 
INDEX X(I) H (I) = ?0 V(I) G (I) >= ?0 U(I) 

1 0 . 1000000E+01 A0.100000E-03 0.120000E+02 

2 0.1000000E+01 A0.299900E+00 0.000000E+00 

3 0 . 1000000E+01 

CHECKING CONVERGENCE THE NORM OF P= 0 . 924648797898831832E-03 


ITERATION 1: 0 PARAMETER ( 1)= 0 . 10001000E+01 

PAGE 1 OBJECTIVE FUNCTION = 0 . 25499301E+02 FUNCTION EVALUATIONS= 14 





CONSTRAINT EVALUATIONS= 28 

INDEX X (I ) H (I) = ?0 

V(I) 

G (I) >= ?0 

U(I) 

1 

0 . 1000483E+01 


A-.133227E-14 

0 . 119995E+02 

2 

0 . 9992333E+00 


A0.299400E+00 

0 . OOOOOOE+OO 

3 

0 . 1000183E+01 





CHECKING CONVERGENCE THE NORM OF P= 0 . 673599522571675425E-03 


ITERATI0N 2: 0 PARAMETER ( 1)= 0 . 10001000E+01 

PAGE 1 OBJECTIVE FUNCTION = 0 . 25499300E+02 FUNCTION EVALUATIONS= 21 





CONSTRAINT EVALUATIONS= 30 

INDEX X(I) H (I ) = ?0 

V(I) 

G (I) >= ?0 

U (I) 

1 

0 . 1000208E+01 


A- . 222045E-14 

0 . 119995E+02 

2 

0 . 9997833E+00 


A0 . 299400E+00 

0. OOOOOOE+OO 

3 

0 . 9999083E+00 





★ ★★★★★★★★★★★★★★★★★★★★★★★★★★ LEAVING RQOPT *********************** 
THE HESSIAN APPROXIMATION UPPER TRIANGLE 
ROW 1 0 . 500000E+01 0.100000E+01 0.100000E+01 
ROW 2 0.500000E+01 0.100000E+01 

ROW 3 0 . 500000E+01 

PERTURBED VALUE OF PARM ( 1 )= 0.999899999999999997 

ENTER THE NUMBER OF ITERATIONS FOR RQOPT? 

1 


Figure A2.7(d) A Sample of RQSEN for Test Problem 2. 
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************************** entering rqopt ********************** 


ENTERING DRIVER ROUTINE 

ITERATION 0: 0 PARAMETER ( 1)= 0 . 99990000E+00 

PAGE 1 OBJECTIVE FUNCTION = 0 . 25500700E+02 FUNCTION EVALUATIONS= 22 

CONSTRAINT EVALUATIONS= 44 
INDEX X(I) H (I) = ?0 V(I) G (I) >= ?0 U(I) 

1 0 . 9997917E+00 AO . 416638E-07 0.120005E+02 

2 0 . 1000217E+01 AO. 300600E+00 0.000000E+00 

3 0 . 1000092E+01 

CHECKING CONVERGENCE THE NORM OF P= 0 . 000000000000000000E+00 
CONVERGENCE ACHIEVED 

*************************** leaving rqopt ********************** 

CENTRAL DIFFERENCE APPROXIAMATIONS TO 
DF/DP (PART F + U PART G) = -6.99999998889772823 
DF/DP (PART F + DF/DX*DX/DP) = -7.00249983049339675 
DF/DP (CENTRAL DIFERENCE) = -7.00249983051293134 

DX/DP (CENTRAL FINITE DIFFERNECE) 

0 . 20831917E+01 -0 . 21667085E+01 -0 . 91669155E+00 


DV/DP (FINITE DIFFERNECE) 

DU ( 1 ) /DP ( 1) (CENT DIFF) = -0 . 46669746E+01 
DU ( 2) /DP ( 1) (CENT DIFF) = 0 . 00000000E+00 
ACTIVE CONSTRAINT DG ( 2) /DP ( 1) = -0 . 60002999E+01 
PGPP+PXPP*PGPP = DG ( 2) /DP ( 1) = -0. 60002999E+01 

LINEAR ESTIMATE OF WHEN ACTIVE SET WILL CHANGE FOR INCREASE P 
G ( 2) ENTERS THE ACTIVE SET FOR DELTA P = 0.49998E-01 

I.E. WHEN P( 1)= 0 . 10499975E+01 

LINEAR ESTIMATE OF WHEN ACTIVE SET WILL CHANGE FOR DECREASED P 
XMIN ( 2) ENTERS THE ACTIVE SET FOR DELTA P =-0 . 46153E+00 

I.E. WHEN P( 1)= 0 . 53847045E+00 

DO YOU WISH TO CALCULATE THE NEW OPTIMUM FOR A NEW VALUE OF PARM( 1)? 

y 

PERFORMING A PARAMETER STUDY FOR PARM ( 1 ) 

ASSUMING BASE POINT IS STABLE 

BASE POINT VALUE PARM( 1 ) = . 1000000000000000E+01 

DF/DP = -0 . 70000000E+01 

DX/DP ( 1)= 0 . 20831917E+01 -0 . 21667085E+01 -0 . 91669155E+00 

ENTER THE PERTURBATION FOR THE PARAMETER OR 
ENTER 0.0 OR NULL TO EXIT FORM THIS SUBROUTINE? 

? 

.1 

Figure A2.7(e) A Sample of RQSEN for Test Problem 2. 
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THE NEW VALUE OF P ( 1)= 0 . 11000000E+01 

★ ★★★★★★★★★★★★★★★★★★★★★★★★★ ENTERING RQOPT 


ENTERING DRIVER ROUTINE 

ITERATI0N 0: 0 PARAMETER ( 1)= 0 . 11000000E+01 

PAGE 1 OBJECTIVE FUNCTION = 0 . 24893909E+02 FUNCTION EVALUATIONS= 30 

CONSTRAINT EVALUATIONS= 60 
INDEX X(I) H (I) = ?0 V(I) G (I) >= ?0 U(I) 

1 0 . 1208319E+01 A0.208111E-01 0.115333E+02 

2 0 . 7833292E+00 A-. 300030E+00 0.000000E+00 

3 0.9083308E+00 

CHECKING CONVERGENCE THE NORM OF P= 0.219176271755135058 

ITERATI0N 1; o PARAMETER ( 1)= 0 . 11000000E+01 

PAGE 1 OBJECTIVE FUNCTION = 0 . 24742404E+02 FUNCTION EVALUATIONS= 37 

CONSTRAINT EVALUATIONS= 62 
INDEX X(I) H (I) = ?0 V(I) G (I) >= ?0 U(I) 

1 0 . 1045910E+01 AO . 415223E-13 0.104876E+02 

2 0 . 7944067E+00 AO . 235367E-13 0.542741E+00 

3 0 . 1055092E+01 

CHECKING CONVERGENCE THE NORM OF P= 0 . 000000000000000000E+00 
CONVERGENCE ACHIEVED 

★ ★★★★★★★★★★★★★★★★★★★★★★★★A* LEAVING RQOPT ********■******★★★★★★★★ 

ENTER THE PERTURBATION FOR THE PARAMETER OR 
ENTER 0.0 OR NULL TO EXIT FORM THIS SUBROUTINE? 

2 

0.0 

DO YOU WANT TO STUDY FINITE PERTURBATIONS 

ENTER PARAMETER NUMBER OR (-1 OR CTRL Z) TO CALCULATE GRADS? 

2 

-1 

DO YOU WISH TO FIND PARTIALS OF THE DESIGN VARIABLES AND 
LAGRANGE MULTIPLIERS W.R.T. PARAMETER (NUMBER OR -1 TO END)? 

2 

-1 

Ready; T=2. 81/4.12 11:23:12 $2.46 

Figure A2.7(f) A Sample of RQSEN for Test Problem 2. 
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